Gravitational Wave Astronomy: 
in Anticipation of First Sources to be Detected 



The first generation of long-baseline laser interferometric detectors of gravitational waves 
will start collecting data in 2001-2003. We carefully analyse their planned performance and 
compare it with the expected strengths of astrophysical sources. The scientific importance 
of the anticipated discovery of various gravitatinal wave signals and the reliability of theo- 
retical predictions are taken into account in our analysis. We try to be conservative both 
in evaluating the theoretical uncertainties about a source and the prospects of its detection. 
After having considered many possible sources, we place our emphasis on (1) inspiraling 
binaries consisting of stellar mass black holes and (2) relic gravitational waves. We draw 
the conclusion that inspiraling binary black holes are likely to be detected first by the ini- 
tial ground-based interferometers. We estimate that the initial interferometers will see 2-3 
events per year from black hole binaries with component masses 10-15 M©, with a signal- 
to-noise ratio of around 2-3, in each of a network of detectors consisting of GEO, VIRGO 
and the two LIGOs. It appears that other possible sources, including coalescing neutron 
stars, are unlikely to be detected by the initial instruments. We also argue that relic gravita- 
tional waves may be discovered by the space-based interferometers in the frequency interval 
2 X 10"'^ Hz-10~^ Hz, at the signal-to-noise ratio level around 3. 
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1 Introduction 



The goal of this review article is quite ambitious. We want to foretell the first gravitational 
wave signals that will be seen by sensitive detectors, several of which are currently in the final 
stage of construction. The detectors will start collecting data in a couple of years from now. 
Obviously, we present a subjective point of view. It is based on our evaluation of what we consider 
the best theoretical knowledge available today, in conjunction with the expected sensitivity of 
the instruments. Possibly, other authors would regard other sources more promising, and would 
place their bet on something else. It is also possible that our view is biased, because it is partially 
guided by the work that we personally were involved in. We will not be very disappointed if we 
are proved wrong. Nature may have many surprises in store for us. It is important, however, that 
for the first time in the long history of gravitational wave research, the conservative astrophysical 
estimates overlap with the detecting capabilities of real instruments. It is an appropriate time 
to prepare strategies for the search and analysis of signals that appear to be more probable than 
others. 

The general theory of gravitational radiation is well understood and is described in textbooks 
P, |2|, ^. The status of the gravitational wave astronomy has been regularly reviewed |Q, ||, 
including papers in Uspekhi 0, ^, IC]. Here, we will only remind the reader that the 



gravitational waves are an inescapable consequence of Einstein's general relativity and, indeed, 
of any gravitational theory which respects special relativity. Gravitational waves are similar to 
electromagnetic waves in several aspects. They propagate with the velocity of light c, have two 
independent transverse polarisation states, and in their action on masses have analogs of electric 
and magnetic components. Gravitational waves carry away from the radiating system its energy, 
angular momentum, and linear momentum. The gravitational-wave field is dimensionless, and 
its strength is qualitatively characterized by a single quantity — the gravitational wave amplitude 
h. The amplitude falls off in course of propagation from a localized source, in proportion to 
the inverse power of the traveled distance: h oc 1/r. The difficulty of direct detection of 
gravitational waves can be seen from the fact that the expected amplitude h on Earth from 
realistic astronomical sources is exceedingly small, of the order or smaller than 10^^^. The 
conceivable amplitudes from laboratory sources are even smaller than that. This small number 
h enters any possible scheme of detection of gravitational waves and makes the detection difficult 
to achieve. For instance, gravitational waves cause a tiny variation Al of the distance / between 
two free masses: Al = hi. In an interferometer with a 1 km arm- length the variation of 
the distance between the two end-mirrors would be of the order Al = 10~^^ cm. This tiny 
variation is supposed to be measured and distinguished against background noise. However, 
in the cosmos, gravitational waves are an important factor of cosmic evolution. Gravitational 
waves are routinely taken into account in the study of orbital evolution of close pairs of compact 
stars 1 11 ] . The measured secular change of orbital parameters in the binary system of neutron 
stars, which includes the pulsar PSR 1913+16, agrees with the gravitational wave prediction of 
general relativity to within 1% accuracy |12|. For the study of pulsars and this discovery, Hulse 
and Taylor were awarded a Nobel prize in 1993. 

Like any other observational science, gravitational wave astronomy operates with sources, 
detectors, data analysis, and interpretation. In what follows, we devote some discussion to each 
of these notions. However, we are not aimed at reviewing all interesting astrophysical theories 
and all possible signals and detection techniques. We concentrate on sources, which, we believe, 
rest on the most solid theoretical foundation, are scientifically important, and involve minimal 
number of additional hypotheses. To be interesting from the point of view of its detection, 
the source should be sufficiently powerful, should fall in the frequency band of the detector, 
and occur reasonably often during the life-time of the instrument. The frequency range of the 
discussed signals is determined by the frequency intervals of the detectors's sensitivity. The 
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currently operating bar detectors are sensitive at frequencies around 10^ Hz. The ground-based 
laser interferometers are sensitive in the interval 10 Hz - Hz. The space-based laser antennas 
will be sensitive in the interval 10~^ Hz - 1 Hz. The great expectations are related with the 
forthcoming sensitive instruments. The Japanese scientists have already built a 300 m laser 
intergerometer called TAMA. The British-German collaboration is in the phase of completion 
of a 600 m laser interferometer called GEO600 |13]. The French- Italian collaboration is building 
a 3 km interferometer called VIRGO ||l3|.The American project LIGO is building two 4 km 



arm- length interferometers [15|. It is expected that these instruments will become operational 
in 1-2 years. The proposal to build a Laser Interferometer Space Antenna (LISA) [16| has been 
tentatively approved by the European Space Agency and NASA, and the launch may occur 
around the year 2010. There exists also plans for advanced ground-based interferometers, such 
as LIGO-II 0. 

The ability of a given instrument to detect a signal depends on the nature of the signal. 
The burst sources, which accompany cosmic catastrophes, emit gravitational radiation at some 
characteristic frequency during just a few cycles. They have tendency to be inherently powerful, 
but their event rate is very low. It is very unlikely to expect such an event to happen in our own 
Galaxy during, say, a 1-year observational run. To see a few events per year, one needs to survey 
a large (cosmological) volume of space and, hence, to possess a sensitive instrument capable of 
detecting the sources from the edges of this volume. The quasi-periodic astrophysical sources 
are expected to be more frequent than the burst sources, but they produce much weaker signals 
in terms of h. However, the amount of the radiated energy during some long time T may be not 
much smaller than that of a burst source. If one knows, or can model, the temporal structure 
of the signal, one can monitor the detector's output during many cycles that are covered by the 
observation time T. This can make a weak periodic signal not much more diffucult to detect 
than a burst signal. Some rare but reliable astrophysical sources, such as binary neutron stars 
and black holes at their latest stage of evolution, exhibit a kind of quasi-periodic gravitational 
wave signal at the inspiral phase, and more like a burst signal in the last moments of their 
coalescence and merging. The stochastic backgrounds of gravitational waves are typically weak 
and difficult to distinguish from the instrumental noise. However, if one can cross-correlate 
the outputs of two or more instruments, and can do this during a long integration time, the 
stochastic background can also be measured. The fundamentally important relic gravitational 
waves form a sort of a stochastic background. They are the only direct probe of the evolution of 
the very early Universe, up to the limits of the Planck era and Big Bang. It would be extremely 
valuable, even if difficult, to detect relic gravitational waves. 

The balance between the expected scientific payoff and theoretical likelihood of various as- 
trophysical sources versus their detectability by the forthcoming and planned instruments is the 
major thrust of this paper. After having analysed many possible sources of gravitational waves, 
and taking all the factors into account, we place our emphasis on compact binaries (neutron 
stars and black holes) and relic gravitational waves. In fact, we argue that inspiraling black 
holes, formed as a result of stellar evolution, are the most likely sources to be seen first by the 
forthcoming sensitive instruments. Also, we think that relic gravitational waves are likely to be 
detected by the advanced ground-based and space-based laser interferometers. To justify our 
point we go into a great detail in describing compact binary stars and relic gravitons. 

Section |^ is devoted to the formation and evolution of binary systems. Binary stars are as 
numerous as single stars. Binaries emit gravitational radiation at twice their orbital frequency. 
To radiate gravitational waves with large intensity and at frequencies accessible to ground-based 
interferometers, the objects forming a pair should be massive and should orbit each other at 
very small separations — a few hundred kilometers. According to the existing views, these 
massive objects can only be the end-products of stellar evolution — neutron stars and black 
holes. Because of the loss of the angular momentum due to gravitational waves, these binary 
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objects are in the late thousands of cycles at their inspiral phase. They are only tens of minutes 
away from the final coalescence and merging, or, possibly, from another spectacular event, a 
gamma-ray burst. The central question is how many such close systems exist in our Galaxy and 
at cosmological distances. This determines the event rate — the number of coalescence events 
that can occur in a given volume of space during, say, 1 year. A detector, sensitive enough to see 
the most distant objects in this volume, will detect all of them. A detector of lower sensitivity is 
capable of seeing the coalescing systems at shorter distances and, hence, will register a smaller 
number of such systems, or will not be expected to see them at all during a 1-year interval of 
observation. 

In sub-section 2T we review all the observational data on binary neutron stars. Even these 
data alone, allow one to derive some estimates on the rate of neutron star coalescences. So far, 
there is no observational evidence of binaries consisting of a neutron star and a black hole or 
two black holes. However, we certainly do not see all the products of stellar evolution in binary 
systems. We need to take into account the predicitons of a theory which successfuly explains 
the formation and relative abundance of various populations of observed binaries consisting of 
normal stars and neutron stars. Such a theory predicts the existence of close binaries involving 
neutron stars and black holes, as the outcomes of processes along certain channels of the binary 
evolution. It is these channels of evolution that are most important for gravitational wave 
astronomy. 

The sub-section 2.2 is devoted to the population synthesis method of describing the continu- 
ing birth and future fate of binary stars. The purpose of this analysis is to find the statistically 
expected number of massive and sufficiently close binaries, which could be in their final stage of 
inspiral at the present cosmological time. This means that we are interested only in those bi- 
naries whose expected total life-time, from formation to coalescence, is shorter than the Hubble 
time. As usual, the results of evolution depend on initial conditions and on physical processes 
along the evolutionary path. We combine the well-established observational facts with reason- 
able theoretical assumptions. Two parameters are especially important — the kick velocity w 
imparted to a newly born neutron star during a supernova explosion, and the fraction kBH 
of a pre-collapse massive star that goes into a resulting black hole. The formation of a black 
hole can also be accompanied by the impartation of some kick velocity. A large kick velocity 
can either disrupt a binary system — a would-be powerful source of gravitational waves, or, on 
the contrary, to make the binary orbit more eccentric, thus increasing the gravitational wave 
luminosity. In our evolutionary calculations we vary w and ksH in the observationally allowed 
limits. We also take into account the stellar wind and the loss of mass as factors of binary 
evolution. The kick velocity is so important a factor of binary evolution that we devote to its 
analysis a separate sub-section \1.3\ 

The results of the population synthesis are summarised in Section ^. These results are at the 
same time our predictions for the detection rate of various compact binary inspiral signals. In a 
given cosmological volume of space, the estimated event rate for coalescing black holes is about 
10 times lower than that for coalescing neutron stars and neutron star - black hole systems. 
However, since the masses of black holes are significantly larger than the masses of neutron 
stars, they are more luminous gravitational wave sources than pairs of neutron stars. Hence, a 
given detector can observe inspiralling black holes at greater distances than pairs of inspiralling 
neutron stars. We conclude that a network of initial laser interferometers are likely to see black 
hole inspirals more often than the neutron star inspirals, and as often as 2-3 events per year. 

Section ^ is devoted to transient and periodic sources. They include supernovae explosions, 
various unstable modes in rapidly rotating neutron stars, and quasi-normal modes of black hole 
perturbations. All these sources are interesting and potentially detectable. However, we do not 
place them at the beginning of our priority list. The asymmetric supernovae explosions, as well 
as the merging event of binaries can produce powerful bursts of gravitational radiation, but 
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the estimates of their performance rely on factors which are not well understood theoretically 
and do not have much of observational evidence. The merging event will be probably seen as 
a confirming signature of the inspiral phase, but one cannot rely on this event alone. However, 
if we err in our priorities, a special kind of hypernovae explosions can top the list. As for the 
unstable modes in rotating neutron stars, they require quite sophisticated mechanisms of their 
excitation and can be hampered by viscosity and other physical processes. The collision of black 
holes and quasi-normal modes of newly born black holes is an intriguing possibility, but should 
probably be treated as something to be discovered by gravity wave observations, rather than 
reliably calculated on purely theoretical grounds. 

In Section ^ we review stochastic gravitational wave backgrounds of astrophysical origin. 
These are the overlaping signals from many individual sources. The populations of sources that 
we consider include unresolved binary white dwarfs in our Galaxy and at cosmological distances, 
and the population of rotating neutron stars. The detection of these backgrounds would carry 
some scientific information on its own, but it is also necessary to study these sources for another 
reason. These stochastic backgrounds set a confusion limit for detection of more interesting 
signals by space-based and ground-based interferometers. We conclude that the LISA will be 
free of the gravitational wave noise from unresolved binaries at frequencies near and higher 
than 2 x IQ-^ Hz. This noise is mostly from unresolved binaries in our Galaxy, the extra- 
Galactic binaries contribute only about 10% to this noise. Any detected stochastic background 
at frequencies above 2 x 10~^ Hz in the LISA window of sensitivity is expected to be of primordial 
origin. The population of non-axisymmetric rotating neutron stars could potentially blur the 
view of the ground-based interferometers. However, we find that this background is below 
the instrumental noise of initial interferometers, and can possibly present a problem only for 
the advanced LIGO. We estimate that other stochastic backgrounds of astrophysical origin are 
weaker than those that we have considered. 

Section |6| is devoted to relic gravitational waves. In contrast to all other sources, which are 
based on classical physics, the generation mechanism of relic gravitons includes some elements 
of quantum physics. It is the inevitable zero-point quantum oscillations of gravitational waves 
amplified by the strong, variable gravitational field of the early Universe that ends up in a 
stochastic background of relic gravitational waves measurable today. Despite the fact that the 
existence of this gravitational wave signal involves an extra element — quantum physics, it is 
not less reliable than many other sources. The generation of relic gravitons relies essentially only 
on the validity of general relativity and basic principles of quantum field theory. Since the same 
mechanism is thought to be responsible for the generation of primordial density perturbation 
seeding the formation of galaxies, we present a qualitative picture of this mechanism in sub- 



section 6.1, The calculation of the expected relic gravitational wave background is given for a 
class of cosmological models supported by other observations. In particular, we use the data 
on the measured microwave background anisotropies. The results of this analysis are presented 
in sub-section 3^. We find that in the most favorable case, the detection of relic gravitational 
waves can be achieved by the cross-correlation of outputs in the initial ground-based laser 
interferometers. In the more realistic case, the sensitivity of the advanced ground-based and 
space-based instruments will be needed. We also discuss a specific statistical signature of relic 
gravitons, associated with the phenomenon of squeezing. This phenomenon is also known in 
formal quantum mechanics and quantum optics. The signature of squeezing could potentially 
help in further improving the signal-to-noise ratio (SNR). 

The problems of detectability are systematically referred to throughout the paper. However 
a rigorous discussion of detectors and data analysis is concentrated in Section |7| and Section ^. 
Whenever we qualify a source as detectable or undetectable, we base our conclusions on the 
more detailed treatment of these Sections. Section ^ gives a general description of detectors 
and their sensitivity curves. The important notion in the detection of a signal with a known 
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or suspected temporal structure, is the notion of a template. A template allows one to use 
the mathched filtering technique (sub-section |8.1| ) in order to increase the SNR. This method 
will be indispensable in the search for signals from inspiraling binaries. The related issues 
are the practically accesible number of templates, their overlap in the parameter space, the 
computational cost, etc. These issues are important not only for a confident detection of a 
signal, but also for the extraction of astrophysical information from a signal (sub-sections 8.4 - 
8.5) — the ultimate purpose of the gravitational wave astronomy. 

Some mathematical details on the Keplerian motion of a binary system and gravitational 



reaction force are described in Appendix Al. Appendix [A2| contains some technical issues of 
the mass transfer modes and mass loss in binary stars. Appendix gives post-Newtonian 
expressions for energy and gravitational wave flux. The main conclusions of the review are 
formulated in the Abstract. 



2 Astrophysical sources. Close binary neutron stars and black 
holes. 

In this Section we discuss observational and theoretical estimates for the coalescence rate of 
close binary neutron stars and black holes. We start from a review of observational limits on the 
coalescence rate of binary neutron stars. Then, we describe the basics of the population synthesis 
of binary evolution which allows one to predict theoretically the event rates for systems involving 
neutron stars and black holes. The role of the kick velocity in the binary evolution is discussed 



in subsection 2.2. The expected detection rates in the forthcoming sensitive gravitational wave 



detectors are summarised in Section y. 

2.1 Observational limits on the binary neutron star coalescence rate 

What do we know about compact binary stars and the rate of their mergings on observational 
grounds? More than a thousand of single neutron stars (NS) are currently (|6[)observed as radio - 



pulsars (see [18|; new data are being continuously added at http:/ /puppsr.princeton.edu] ). 
In addition, about 30 NS are seen as X-ray pulsars and yet more 100 NS are seen as burst and 
transient X-ray sources. These NS enter binary systems with non-degenerate companions, that 
is, their companions are normal stars rather than neutron stars or black holes. Only six NS 
are known to enter binary systems with another NS as a secondary component 0. All these six 
systems belong to binary radio-pulsars. The systems and some of their parameters are listed in 
Table |l[ Orbital periods are given in days and masses in units of the solar mass Mq . Three of 
these systems (namely, B1913+16, B1534+12 and B2127+llc) are close enough to merge due 
to GW emission in a time interval shorter than the Hubble time tn- We loosely refer to binaries 
as coalescing or merging binaries if their expected life-time up to coalescence, t^oah is shorter 
than tn- For numerical estimates we use the value ~ 12Gyr = 12 x lO^yr. 

Much less is known about black holes (BH). A dozen of BH candidates participate in binary 
systems with non-degenerate companions. They are observed as persistent X-ray sources (like 



Cyg X-1) or X-ray transients (mostly X-ray Novae) (see |21] for a review). Neither single BH nor 
BH forming a binary with radio-pulsar or another BH have been found so far. Parameters of the 
BH candidates in binary systems are listed in Table Note that according to these data, the 
mean BH mass is Mbh — S.SMq, i.e. notably higher than a typical NS mass Mns — 1.4Mq. (Of 
course, we mean a black hole in astrophysical sense, i.e. as a highly compact gravitating object 
of certain mass. The presence or absence of the event horizon is irrelevant for our discussion.) 



For a recent summary of NS mass determination see [22|. 



^New binary pulsars are found in recent pulsar surveys (see e.g ^). However, a reliable determination of 
the component masses is only possible after sufficiently long-term observations. 
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Table 1: Binary PSR with NS secondaries (data from [p3[) 



PSR 


Pb{d) 


e 


M1 + M2 


Ml 


M2 


^coalesi 


J1518+4904 


8.634 


0.249 


2.62 






>3.610i2 


B1913+16^ 


0.323 


0.617 


2.8284 


1.44 


1.39 


1.010*^ 


B1534+121 


0.420 


0.274 


2.6784 


1.34 


1.34 


1.0 10^ 


B2127+llc^'2 


0.335 


0.681 


2.712 


1.35 


1.36 


8.010^ 


B2303+46 


12.340 


0.658 


2.60 






>1.610^2 


B1820-113 


357.762 


0.795 








>2.410^5 



Coalescing binary pulsars 
^ Binary pulsar in a globular cluster 
^ The secondary companion may not be a NS 



Table 2: BH Candidates [|T| 



System 



Spectral class 


Porb, d 


fv{m),MQ 


irix, 




my,MQ 


09,7 lab 


5.6 


0.23 


7- 


18 


20-30 


B(3 6)II-III 


1.7 


2.3 


7- 


11 


3-6 


0(7-9)111 


4.2 


0.14 


4- 


10 


18-25 


K(5-7)V 


0.3 


3.1 


5- 


17 


~0.7 


KOIV 


6.5 


6.3 


10- 


-15 


0.5-1.0 


K(3-5)V 


0.4 


3.01 


9- 


16 


0.7-0.8 


K(3-7)V 


0.3 


5.0 


5.3- 


-8.2 


~0.7 


M(0-4)V 


0.2 


0.9 


2.5- 


-5.0 


~0.4 


F5IV 


2.6 


3.2 


4- 


-6 


~2.3 


K3 


0.7 


4.0 


5- 


-7 


~0.8 



Cyg X-1 
LMC X-3 
LMC X-1 
A0620-00 
GS2023-F338 
GSR1121-68 
GS2000-h25 
GRO J0422-h32 
GRO J1655-40 
XN Oph 1977 

Cyg X-3 ? 



Mean Value of the BH mass 



■'8.5Mp 







There are two types of estimates of the binary NS coalescence rate. The estimates of the 
first type are derived directly from observations (see Table while the estimates of the second 
type are inferred from theory of binary stellar evolution (see Table 0) . We will consider each of 
these estimates in turn. 

The estimates of the first type are based on the data on three binary radio-pulsars which 
should merge within the Hubble time (Table |l|). These estimates use the following argumen- 
tation. The average coalescence time for these pulsars is approximately 3 x 10^ years. So the 
binary NS merging rate based on these 3 pulsars would be approximately once per 100 million 
years. As we observe only about 1% of the Galactic volume, a lower limit for the binary NS 
merging rate becomes one every million years [24|. In fact, this estimate was formulated at the 
time when only two of the presently known three coalescing binary radio-pulsars were known. 
Taking into account the spatial distribution of pulsars inside the Galaxy and the fact that a 
typical radiopulsar switches-off long before the coalescence, the lower limit for NS merging rate 
can be increased by almost an order of magnitude [^] , thus reaching 10~^ per year. 

An interesting upper limit, the so-called "Bailes limit", was derived from independent argu- 
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Table 3: Observational estimates of binary NS coalescence rate. 



Author (s) 



Coalescence rate (yr ^] 
1/10^ 
1/10^ 
3/10^ 
8/10^ 

< 1/10^ 

< 1/10^ 



Phinney 1991 
Narayan et al 1991 [^] 
Curran, Lorimer 1995 |26] 
van den Heuvel, Lorimer 1996 



"Bailes limit" 1996 |28| 
Arzoumanian et al. 1999 



ments [p8|| . It was noted that the properties of the pulsars in the three coalescing binary radio 
pulsars (most of all, their surface magnetic fields) are quite different from those found in ordi- 
nary single radio pulsars. Since the number of single radio pulsars is about 1000, it is estimated 
that the radio pulsars similar to those residing in merging binary NS should be formed al least 
~ 1000 times rarer than single radio pulsars. Taking the birth rate of single pulsars from a large 
sample of known pulsars, Bailes proposed to put an upper bound on the birth rate of binary NS 
as the (formation rate of single pulsars) x (number of pulsars with ordinary properties among 
binary radio-pulsars) = (1/60 yr)x (1/1000)~ 2 x 10~^yr~^. 

One should note, however, that in both the estimates - the one based on statistics of coalesc- 
ing binary radio pulsars, and the other based on the Bailes' limit - suffer from selection effects. 
They depend on the pulsar distances (in some cases known not better than up to a factor of 2) 
0, on the characteristic pulsar life-time (known not better than up to an order of magnitude), 
and on the differences in the properties of single and binary pulsars. So, one cannot infer from 
observations an absolutely reliable estimate of the binary NS merging rate. Indeed, a recent 
re-assessment of the Bailes limit [^] taking into account the current pulsar numbers and the 
reduction in search sensitivity to short orbital period binaries gave 10~^ yr~^ for the upper limit 
of binary NS Galactic merging rate. An alternative way of deriving the upper limit based on 
empirical pulsar birth rate and theoretical understanding of binary NS formation was used in 
|l3l| to yield a few mergers per 10^ years. The cited papers clearly demonstrate that (1) there 
is a steady tendency to increase the empirical upper limit of the binary NS coalescence rate 
and (2) various selection effects generic to radio pulsar surveys and lack of detailed knowledge 
of Galactic pulsar population properties still prevent us from the derivation of a fully reliable 
estimate. 

2.2 Population synthesis of coalescing binary NS and BH 

Now we turn to estimates partially based on theoretical grounds. The merging rates of binary 
compact stars have been calculated by different independent research groups, mostly with the 
help of the population synthesis numerical simulations (see Table Q). The reliability of these 
results depend on whether the binary evolution scenarios properly reflect different aspects of the 
real, observed populations. It is encouraging that these independent calculations yield similar 
results. 

Theoretical estimates of double NS coalescence rate are systematically higher than the ob- 
servational ones by, on average, an order of magnitude. This does not mean that the estimates 
are in conflict with each other. The main reason for the discrepancy is that the observational 
estimates directly refer only to the merging rates of binary NS, in which one of the components 
is a radio pulsar. The participation of a pulsar in a binary system is of course not a necessary 

^For example, recent observations of PSR 1534+12 |30j suggest a distance which is two times larger than was 
previously thought, so the "observational estimate" of binary NS mergings should be decreased by a factor > 2. 
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condition for the system to be interesting from the point of view of gravitational wave astronomy. 
For example, the neutron stars could be born with weak magnetic fields and/or slow rotational 
periods, and, hence, they would never manifest themselves as radio pulsars. Theoretical cal- 
culations provide a broader range of estimates because they depend on several evolutionary 
parameters which are not well known. However, the population synthesis calculations are the 
only way to estimate the coalescence rates of, so far unobserved, compact binaries consisting 
of two BH or a NS and a BH. These systems are of a great importance for gravitational wave 
astronomy. We will consider below the population synthesis for all possible pairs: NS+NS, 
NS+BH, and BH+BH. 



2.2.1 Basics of population synthesis 

Coalescence rates of different types of binary compact stars are calculated using modern theory 



of binary star evolution (e.g. |4C] and references therein). A full description of the method 



can also be found in [38|. Some key formulas for binary system evolution are summarized in 



Appendix A2 



Binary stars are formed with different initial masses, semi-major axes, eccentricities, etc.. 
These initial parameters are drawn from certain distribution laws. Also, there are some other 
physical parameters important for binary evolution, such as the efficiency acE of the angular 
momentum removal in the common envelope stage, or kick velocity distribution /(w) for newly 
born neutron stars. This means that in order to calculate the expected rate of binary star 
mergers, we need to derive the number of binaries formed in all appropriate regions of the 
parameter space, and then integrate over these regions. Some distributions can be more or 
less accurately deduced from astronomical observations. To these belong the initial stellar 
mass function and the distribution over binary semi-major axis. Distributions of other physical 
parameters are being adopted on theoretical grounds. We perform evolutionary calculations 
using the "Scenario Machine" code - a version of the Monte-Carlo method (see [^8| for review). 
In a typical numerical experiment, some ~10^ binary evolutionary tracks with different initial 
conditions are calculated. Similar approaches in the binary evolution studies are being used by 
other groups (e.g. p^]) and are commonly named "population synthesis" methods. 



2.2.2 Initial binary parameters 

The initial components of a binary system are taken as zero age main-sequence stars. The initial 
parameters which determine the subsequent binary evolution are: the mass of the primary 
component Mi, the binary mass ratio q = M2/M1 < 1, and the orbital separation a. For 
sufficiently close binaries, which are capable of producing a merging NS or BH, we assume that 



Table 4: Theoretical estimates of binary NS coalescence rate. 



Author (s) 



Coalescence rate (yr 
TJW- 



-i/i(F 

1/10^ 
1/10^ 
S/IOM/IO"^ 
< 3/10^ 
3/10^ 

' 1/10^-3/10^ 



Clark et al 1979 
Lipunov et al 1987 



Hills et al 1990 |4| 
Tutukov, Yungelson 1993 |35] 
Lipunov et al 1995 



Portegies Zwart, Spreeuw 1996 |37] 
Lipunov et al 1996 || 
Portegies Zwart, Yungelson 1998 Il39| 
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the orbits are circular. This assumption is justifiable since the tidal interaction between the 
components is effective enough to circularize the orbit. 

The distribution of binaries over initial orbital separations is partially known from observa- 
tions 

/(logo) = const, for a such that 
max {10 Rq, Roche Lobe (Mi)} < a < 10^ Rq. 



7. (1) 



We assume that the distribution of masses of the primary (more massive) components obeys 
the Salpeter mass function found for the birth rates of main-sequence stars in the solar vicin- 

ity m-. 



f 



Ml \ / M_ 



-2.35 



The observed star formation rate in the Galactic disk relates to this distribution as 



0.9 yr~V , 0.1 Mq < Mi < 120 , 



dN 

dMidt ' VMe 



(3) 



Assuming that 50% of the total number of stars in the Galaxy reside in binaries, this distribution 
law predicts one massive star (Mi > IOMq to be able to produce a compact remnant) to form 
in a binary system, approximately every 60 years. This estimate agrees with the binary birth 
rates derived from observations of binary stars |43|. 



The initial mass ratio g in a binary is crucial for its subsequent evolution |44] because it 
determines the mode of the first mass transfer between the components. The initial distribution 
over q has not been reliably derived from observations due to various selection effects. Usually, 
one makes the 'zero-order assumption', according to which the mass ratio distribution is flat, 
i.e. low mass ratio binaries are formed as frequently as those with equal masses (e.g. [p0|]): 

dN 

— = const , q = M2/M1 < 1 . (4) 
dq 

In the calculations presented below we used this prescription. The effect of the initial q distri- 
bution on binary pulsar statistics was studied in ref. fl^ . 

2.2.3 Neutron star kick velocities 

In the course of evolution of massive binary stars, one or two neutron stars can form. One of 
the most important parameters affecting the eventual binary NS coalescence rate is the kick 
velocity w imparted to a NS at its birth. There exists a plenty of observational evidence for 
the kick velocities. The impact of a kick velocity ~ 100 km/s explains the precessing binary 



pulsar orbit in PSR J0045-7319 |46]. The evidence of the kick velocity is seen in the inclined. 



with respect to the orbital plane, circumstellar disk around the Be star SS 2883 — an optical 



component to the binary pulsar PSR B1259-63 [47|. One more direct evidence comes from 
observations of the geodetic precession in the famous binary pulsar PSR 1913+16 |48, 
All these observations indicate that in order to produce the observed mis-alignment between 
the orbital angular momentum and the neutron star spin, a component of the kick velocity 
perpendicular to the orbital plane is required. A non-zero kick velocity is also required in order 
to properly explain the observed pulsar velocity distribution 0. 

Most likely, the origin of the kick velocity should be attributed to the asymmetry in the 
supernova explosion. Astrophysical evidence for the existence of a substantial kick velocity 



during supernova explosions was discussed in |52] and recently summarized in However, 



the nature of concrete physical mechanisms giving rise to the kick velocity is still unclear (see 
''For an alternative point of view see ^] and for its criticism ]5l| . 
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1 55 1 for recent review). A promising possibility capable of producing small and moderate kick 
velocities (up to 100 km/s) involves asymmetric neutrino emission in strong magnetic field of a 
newly born NS [||, 

Spatial velocities of NS are usually derived from direct observations of proper motions of 



single radio pulsars |58, p^, |5^. Or, with more uncertainty, from the observed offsets in the 



positions of young pulsars relative to the centers of their associated supernova remnants (e.g. 
pO| ). Both methods determine, however, only the component of the space velocity which is 
transversal to the line of sight. One should also bear in mind the existing uncertainty in pulsar 
distances, which affects the evaluation of the kick velocity. On average, the uncertainty in the 
distance scale is ~ 30% [P], and can be as much as a factor 2 in individual cases. In general, 
the observed distribution of transversal pulsar velocities is recognized to have, (1) a high mean 
value (~200-350 km/s) and (2) a broad shape with a high-velocity tail up to 1500 km/s. 

It is a difficult problem to derive the intrinsic kick velocity distribution from these data. If 
all the pulsars, presently seen as single pulsars, took their origin from single massive stars, their 
velocity distribution would have exactly reflected the initial kick velocity distribution. This is 
because single massive stars have very small (of order 10 km/s) spatial velocities. However, 
when a supernova explosion occurs in a binary system, which gets disrupted as a result of the 
explosion, the neutron star can acquire a substantial space velocity, equal to the orbital velocity 
of the progenitor, even without any additional kick velocity. If a non-zero kick is present, it 
becomes practically impossible to solve analytically the inverse problem for the intrinsic kick 
velocity distribution from the observed pulsar velocities. 

So, the only way to check the very assumption of the non-zero kick is to find numerically the 
pulsar velocity distributions arising from various theoretical kick distributions, and compare the 
results with observations. This is usually done by the Monte-Carlo simulations of binary evolu- 
tion. Presently, there is no general agreement with regard to the form of the kick distribution. 
A Maxwellian distribution for ii; = |w| 

fm{w) oc w"^ exp(-w^/wo) , (5) 
with wq = 190 km/s, was used by Hansen and Phinney p2| in order to fit the observed pulsar 



velocities. A different form of the kick velocity distribution was suggested in [45|. It was found 
that this latter distribution fits well with the observed 2D pulsar velocity distribution given by 
Lyne and Lorimer [^]. In contrast to the Maxwell-like distribution, the proposed distribution 
function has a power-law shape: 

, n 1^ (|w|/^0)°-^^ 

^^^^'"^'^ (l + (|w|/u;o)6-^^)°-^ 

and assumes wq « 400 km/s. In our calculations described below, we have used this form of 
the NS kick velocity distribution, where wq is treated as a free parameter. A study of the kick 
velocity effects on the binary neutron star merging rate can be found in |6S] and [^]. We give 



a detailed analysis of these effects in subsection 2.3 



2.2.4 Binary neutron star formation and merging 

We are interested in evolutionary tracks which lead to the formation of a pair of coalescing 
NS. Detailed studies of possible evolutionary channels which produce merging binary NS can 
be found in the literature (e.g. |6^, ^ |6^)- Usually, the evolutionary analysis is being 
done in the following order: one starts from the observed parameters of the binary and tries 
to deduce the parameters of the supernova progenitor and then, the initial binary masses and 
orbital separation. In contrast, the Monte-Carlo population synthesis method, which we apply, 
evolves a trial binary and looks for appropriate results by changing the initial parameters within 
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Figure 1: Evolutionary track of a massive binary star leading to the formation and coalescence 
of two NS. 



their distributions. One of such typical tracks calculated by us is sketched in Fig. |l], and we will 
explain it in detail. 

Close NS binaries originate from two sufficiently massive main-sequence stars separated 
by a moderate distance of order 100 solar radii (1-st row in Fig. [l|). To get a NS in the 
course of evolution, the mass of the progenitor star must be larger than wlOM© initially or, 
in any case, taking into account a possible mass transfer in close binary, the mass should be 
~1OM0 during the stage of nuclear burning. The more massive the primary star is, the faster 
it evolves. For main-sequence stars, the time of core hydrogen burning is tnud M~^. The 
star burns out its hydrogen in its central parts, so that a dense central helium core with a mass 
Miif, ~ O.1(M/M0)^'^ forms by the time when the star leaves the main-sequence. The outer 
shell expands and the star moves towards the red super-giant region in the Hertzsprung-Russell 
diagram. At some stage of its evolution the star fills out its Roche lobe [Eq. ( [A42| )] (3-rd row in 
Fig. HI). The hydrogen envelope starts outflowing onto the secondary, less massive star, which 
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still resides on the main-sequence. The primary star is being continuously stripped off of its 
hydrogen envelope, until a naked helium core emerges. This core can be observed as a hot 
compact helium star, or, for more massive stars, as a Wolf-Rayet star with intensive stellar 
wind (5-th row). 

While the mass of the primary star reduces, the mass of the secondary star increases, since 
the mass transfer at this stage is thought to be quasi-conservative. For not too massive main- 
sequence stars, M < 20 Mq, no significant stellar wind mass loss occurs which could, otherwise, 
remove matter from the binary. The secondary star acquires a large angular momentum due to 
the infalling material, so that its outer envelope can be spun up to an angular velocity close to 
the limiting (Kepler orbit) value. Such massive rapidly rotating stars are observed as Be-stars. 
During the conservative stage of mass transfer, the semi-major axis of the orbit first decreases, 
reaches a minimum when the masses of the binary components become equal to each other, 
and then increases. This behavior is dictated by the angular momentum conservation law [Eq. 
( A26| )]. After the completion of the conservative mass transfer, the initially more massive star 
becomes less massive than its initially lighter companion. The parameter q = M2/M1 becomes 
larger than 1. In a short time, typically ~ 10% of the hydrogen burning time, the nuclear 
evolution of the helium star is completed and, provided its mass is larger than 2-3 Mq, it 
explodes as a core-collapse supernova type-II leaving a neutron star as its remnant. 

Even for asymmetric supernova explosions, most of such binaries do not get disrupted. This 
is because the mass ratio q = M2/M1 of the pre-supernova binary becomes generally high, 
q ~ 3-5. After the first SN explosion, the binary system consists of a Be-star and a NS in an 
elliptical orbit (8-th row). Orbital evolution following the SN explosion is described in more 
detail in an Appendix [Eqs. ( |A35| - |A40| )]. 

Be-stars have very rapidly rotating envelopes but in other respects they do not differ from 
ordinary main-sequence stars. After the completion of hydrogen burning in the core, a Be-star 
starts expanding until it fills out its Roche lobe while passing through the periastron of an 
elliptical orbit (10-th row). This initiates the second episode of the mass transfer, which takes 
place on the thermal scale of the Be-star, typically ~ lO~^M0/yr. However, this mass transfer 
is qualitatively different from the first one, since the mass transfer is now on to a compact star. 
Once the accretion rate exceeds the value which provides the luminosity equal to the Eddington 
luminosity limit near the NS surface (~ 10~^MQ/yr), the NS cannot accrete all the infalling 
matter. The so-called common envelope stage arises (11th row) during which the neutron star 
finds itself inside quite dense outer layers of the companion star. Numerical hydrodynamic 
calculations |67, 38 1 show that the dynamical friction of the orbiting NS leads to an efficient 
transfer of the orbital angular momentum to the common envelope, thus dispersing it on a very 
short timescale (typically, 10^ — 10"^ years). The semi-major axis of the binary system reduces 
dramatically [Eq. ( |A41[) ], which results in the formation of a close binary system consisting of 
a NS and a WR star (12th row). Alternatively, the NS can sink into the center of its red giant 
companion (the so-called Thorne-Zytkow object; not shown in this Figure). 

In a short time (^ 10^ years), the companion WR star explodes as a supernova type lb, 
thus producing a second neutron star. During this explosion, the system is more likely to be 
disrupted than during the first SN explosion, since the exploding star is now more massive than 
its companion. Surviving systems form close high-eccentric NS binaries, similar to the NS binary 
PSR 1913-1-16. Orbital parameters of such binaries change exclusively due to the emission of 
gravitational waves (see section Al.4). If the neutron stars are close enough, they coalesce in a 
time shorter than tu- 
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2.2.5 Black hole formation parameters 

So far, we have considered the formation of individual NS and their binaries. It is believed that 
very massive stars end up their evolution with the formation of stellar mass black holes. We 
will discuss now the formation of an individual black hole. 

In the analysis of BH formation, new important parameters appear. The first one is the 
threshold mass Mcr beginning from which a main-sequence star, after the completion of its 
nuclear evolution, can collapse into a BH. This mass is not well known, different authors assume 
different values: van den Heuvel and Habets |69] — 40Mq; Woosley et al. |Q — 60Mq; Portegies 



Zwart, Verbunt, Ergma [^] — more than 2{)Mq. A simple physical argument usually cited in the 
literature is that the mantle of the main-sequence star with M > Mcr ~3OM0 is bound before 
the collapse with the binding energy well above 10^^ ergs (typical supernova energy observed), 
so that the supernova shock is not strong enough to expel the mantle. 

The second parameter is the mass Mbh of the formed BH. There are various studies as for 
what the mass of the BH should be (e.g. |72, 73, 74, |7^). In some papers a typical BH mass 



was found to be not much higher than the upper limit for the NS mass (Oppenheimer-Volkoff 
limit ~1.6-2.5M0, depending on the unknown equation of state for the neutron star matter) 
even if the fallback accretion onto the supernova remnant is allowed [^]. However, observations 
strongly indicate much higher masses of BH candidates, of the order of G-IOM© (see Table ^) . To 
obtain such BH masses, it is sometimes assumed [^] that Mcr ~80Mq. Recently, a continuous 
range of BH masses up to 10-15 Mq was derived in calculations |7^. Since the present day 
calculations are still unable to reproduce self-consistently even the supernova explosion, we have 
parameterized the BH mass Mbh by the fraction of the pre-supernova mass M^, that collapses 
into BH: ksH = Mbh/M^,. In fact, the pre-supernova mass M,, is directly related with Mcr, 
but the form of this relationship is somewhat different in different scenarios for massive star 
evolution. According to our parameterization, the minimal BH mass can be Mj^^ = kBnM^, 
where M* itself depends on Mcr- We have varied kBH in a wide range from 0.1 to 1. 

The third parameter, similar to the case of NS formation, is a possible kick velocity wbh 
attributed to a newly formed BH. In general, one expects that a BH should acquire a smaller 
kick velocity than a NS, as black holes are more massive than neutron stars. In our calculations 
we have adopted the relation 

WBH _ - Mbh _ 1 - kBH 
WNS ~ - Mov ~ 1 - Mov/M^ ' 

where Mqv = 2.5M0 is the maximum NS mass. When Mbh is close to Mov, the ratio 
wbh /wj\fS approaches 1, and the low-mass black holes acquire kick velocities similar to those 
of neutron stars. When Mbh is significantly larger than Mov, the parameter knn = 1, and 
the BH kick velocity becomes vanishingly small |^. As we show below, the allowance for a quite 
moderate wbh strongly increases the coalescence rate of binary BH. A recent analysis of space 
velocities of some BH candidates did not reveal the need for a non-zero wbh [^6|. However, 
other studies show that some kick velocity can arise during the BH formation, and its presence 



does not contradict the observational data |77|. From a theoretical point of view, the presence 



of a moderate kick velocity imparted to a BH during its formation seems very plausible [ ]75[ . 

2.2.6 Binary black hole merging with wbh = 0: typical example 

We begin from the simplest assumption wbh = 0. The more realistic cases wbh 'wiH be 



considered in subsection 2.3.2. In contrast to NS-I-NS binaries, the BH-I-BH and BH-l-NS binary 



*Other possible relationships between wbh /wns have also been checked, but their different forms do not affect 
the results significantly. 
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Figure 2: Evolutionary track of a massive binary star leading to the formation and coalescence 
of two BH. The low stellar wind mass loss scenario is used. 



systems have not been observed so far. There is no way of recovering from observations the range 
of progenitors for such binary systems. We can only apply the population synthesis method and 
derive theoretically the parameters of all the binaries, including the BH+BH and BH+NS pairs, 
that should be produced at the end of the evolution of very massive binary stars. 

In addition to the evolutionary uncertainties existing for stars evolving to binary NS systems, 
new uncertainties arise for very massive stars with initial masses M > 40Mq. First of all, a 
large mass-loss via stellar wind is observed for such stars. According to current views, a massive 
single star can lose more than a half of its initial mass already on the main-sequence. Further 
rapid mass decrease is expected during a helium star stage. There is no general agreement as 
to how to describe the mass loss of a massive star. Yet, one can consider two extreme cases for 
the mass loss via stellar wind: a slow mass-loss and very fast one. Since the exact description 
of the stellar wind mass loss is not known, we have considered both options in our numerical 
simulations (for more details, see ||6^). A typical evolutionary track that leads to the formation 
of coalescing binary BH system, assuming a low mass-loss scenario, is shown in Fig. |2|. 

Before proceeding further with evolutionary calculations, we want to explain qualitatively 
how a merging BH binary can form even in the framework of an extremely high mass loss. The 
high mass loss makes the binary system wider, according to Eq. ( A2g| ), so it should have been 
quite tight before the phase of an active mass-loss. Consider two already existing black holes, 
each of 10 Mq, in a circular orbit. The orbital separation should be < 20Rq for the binary to 
merge within the tff (see Fig. |l^). This means that the radius of a Wolf-Rayet (helium) star 



17 



that collapsed second, in the course of a binary evolution, to form a BH should have been less 
than 10 Rq. The mass of the pre-collapse star could not be smaller than 10 Mq, since ksH < 1- 
Such massive helium stars have very small radii (~ IRq) and do not expand too much before 
the collapse, so the requirement R < IORq is fulfilled. The life-time of a massive helium star is 
about 10^ years, since it loses mass at a high rate of 10~^ Mq/yt. The star can lose a sizable 
fraction, maybe a half, of its mass before the collapse. Thus, we will be dealing with a 20 
Mq helium star in pair with the first 10 Mq BH in a circular orbit of radius a ~ 13i?0. (We 



applied Eq. ( A29 ) to calculate the radius of the resulting orbit.) Note that the 8 Rq Roche lobe 
of a 20 Mq helium star is still quite large. To form such a close WR+BH binary, a common 
envelope stage is needed. The 20 Mq He core corresponds to at least 55 Mq main-sequence 



star, as follows from Eq. ( A44| ). According to the models by Schaller et al. |78], a massive 



star loses about a half of its initial mass on the main-sequence, so to form a common envelope 
with 10 Mq BH the star has to lose ~ 25M0 while on the main-sequence. This means that the 
common envelope stage should have started with a 30 Mq red super-giant filling its Roche lobe 
and having a 10 Mq BH as a companion. The mass ratio in such a system is high enough for 
the common envelope to develop. The orbital separation at the common envelope stage should 
have decreased by 6-12 times, according to Eq. ( A43| ), depending on the parameter acE and the 



exact value of the red giant mass. So, before the common envelope stage, the orbital separation 
should have been ~ ISO-Rq. The orbit should be somewhat smaller than this (i.e., about 120 
Rq or less) when the first BH forms, because of the strong wind from the red giant and loss of 
the total mass. And in order to collapse first, the mass of the primary star must have been at 
least 60 Mq. Assuming isotropic stellar wind and using again Eq. ( [A29D , we conclude that the 



initial system could have widened at most (60 + 55)/(10 + 30) ~ 3 times since the time of its 
formation, i.e. the initial separation of the progenitor binary should be larger than 40 Rq. The 
initial separation of 50 Rq is sufficient enough to harbour two 60 Mq stars since their radii are 
less than 20 Rq on the main-sequence. Even though such initially close massive binaries are 
rare, they should exist. Thus, we see that some fraction of massive binary stars should have 
ended up as sufficiently close pairs of black holes. 

2.3 Effects of the kick velocity 

The picture outlined above changes if a non-zero kick velocity is present in the process of 
formation of a NS or a BH. This, in turn, has a significant effect on the expected rate of 
compact binary mergings, which is of primary interest for our study. In general, the formation 
of a compact object (NS or BH) is accompanied, both, by a mass loss from the system and 
by a kick velocity. The effects of kick velocity during supernova explosions were considered in 



many papers (see e.g. |79, B0| and also A2). The general formulae for the condition of system's 
disruption and for parameters of the resulting elliptical orbit, if the system remains bound, 
are derived in A2, see Eqs. ( A35| ), ( |A36| ) , ( |A38 ) . Here we will present qualitative arguments 



enabling the reader to see the main consequences of a non-zero kick velocity. We restrict our 
attention to circular orbits and assume equal probabilities for all possible orientations of the 
kick velocity vector w. We argue that a moderate (not too large) kick velocity increases the rate 
of binary mergings. This happens because a moderate kick velocity does not change too much 
the likelihood of the system's disruption, but, at the same time, always makes the periastron of 
the resulting elliptical orbit smaller than it would have been without a kick. As a result, some 
of the binaries, whose coalescence time without a kick would be longer than the Hubble time, 
now get a chance to merge in a time shorter than tn- This increases the number of detectable 
gravitational wave sources. 
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2.3.1 Effect of the kick velocity on the disruption of a binary system 

The collapse of a star to a BH, or its explosion leading to the formation of a NS, are normally 
considered as instantaneous. This assumption is well justified in binary systems, since typical 
orbital velocities before the explosion do not exceed a few hundred km/s, while most of the 
mass is expelled with velocities about several thousand km/s. The exploding star Mi leaves the 
remnant Mc, and the binary loses a portion of its mass: AM = Mi — M^. The relative velocity 
of stars before the event is 



V- = ^G(Mi + M2)M. (8) 
Right after the event, the relative velocity is 

Vj = Vi + w. (9) 

Depending on the direction of the kick velocity vector w, the absolute value of Vj varies in the 
interval from the smallest Vf = \ Vi — w\ to the largest Vf = Vi + w. The system gets disrupted 
if Vf satisfies the condition (see A2): 




(10) 

V ^ 

where x = (Mi + M2)/{Mc + M2). 

Let us start from the limiting case when the mass loss is practically zero (AM = 0, x = 1); 
while a non-zero kick velocity can still be present. This is a model for a BH formation with 
ksH = 1- It follows from Eq. (|lO|) that, for relatively small kicks, w < {V^ — l)Vi, the system 
always (independently of the direction of w) remains bound, while for w > (\/2 + 1)^^ the 
system always unbinds. By averaging over equally probable orientations of w with a fixed 
amplitude w, one can show that in the particular case w = Vi the system disrupts or survives 
with equal probabilities. If Vf < Vi, the semi-major axis of the system becomes smaller than 
the original binary separation, af < Ui (see Eq. ( |A35 )). This means that the system becomes 



more bound than before, i.e. it has a greater negative total energy than the original binary. If 
Vi < Vf < ^/2Vi, the system remains bound, but a/ > Oj. For small and moderate kicks w <Vi, 
the probabilities for the system to become more or less bound are approximately equal. 

In general, the binary system loses some fraction of its mass AM. For a BH formation 
this corresponds to ksH < 1- In the absence of the kick velocity, the system remains bound 



if AM < M/2 and gets disrupted if AM > M/2 (see |A2D . Clearly, a properly oriented kick 
velocity (directed against the vector Vj) can keep the system bound, even if it would have been 
disrupted without the kick. And, on the other hand, an unfortunate direction of w can disrupt 
the system, which otherwise would stay bound. 

Consider, first, the case AM < M/2. The parameter x varies in the interval from 1 to 2, 
and the escape velocity Vg varies in the interval from \^Vi to Vi (see |A^ ). It follows from Eq. 



( A39| ) that the binary always remains bound if w < Ve — Vi, and always unbinds if w > Ve + Vi. 



This is a generalization of the formulae derived above for the limiting case AM = 0. Obviously, 
for a given w, the probability for the system to disrupt or become less bound increases when 
AM becomes larger. Now turn to the case AM > M/2. The escape velocity of the compact 
star becomes Ve < Vi. The binary is always disrupted if the kick velocity is too large or too 
small: w > Vi + Ve or w < Vi — Ve. However, for all intermediate values of w, the system can 
remain bound, and sometimes even more bound than before, if the direction of w happened 
to be approximately opposite to Vj. A detailed calculation of probabilities for the binary's 
survival or disruption requires integration over the kick velocity distribution function /(w) (see 

e.g. m- 
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2.3.2 Effect of the kick velocity on coalescence rate of compact binary systems 

Here we consider binary systems that were not disrupted during tlie formation of a compact 
object. The parameters aj and e of the resulting elliptical orbit are defined by Eqs. ( [A35| ), 
( A36| ). The distance of the closest approach between the stars is given by the orbit's periastron 
Gp = af(l — e). It follows from Eqs. ( |A35| ), ( A36| ) that Up = ai in the absence of kick velocity. 



The importance of the kick velocity w ^ lies in the fact that, although the semimajor axis o/ 
can increase or decrease under the action of the kick, the periastron distance always becomes 
smaller: Op < a^. This relationship follows from the combination of Eqs. ( A35 ), ( A36| ) plus 



the requirement that the system remains bound, i.e. the quantities participating in Eq. ( A38| ) 
satisfy the opposite inequality. The decrease of the periastron distance plays an important role 
in the subsequent evolution of the binary, which consists now of a newly born compact star and 
its companion. 

Consider, first, a normal star as the companion. Since the kick has diminished the perias- 
tron distance, as compared with the no-kick case, the normal star, while passing through the 
periastron, will fill out its Roche lobe in a shorter time, than it would do in the absence of the 
kick. After the tidal circularization of the orbit, a tighter binary is formed. Accordingly, the 
subsequent common envelope stage makes the binary tighter than it would otherwise do (see 
Eq. ( A43| )). As a result, the final binary system, consisting of two compact objects, will coalesce 



due to GW radiation in a shorter time (see Eq. (|A22| )). In other words, some of the binaries, 
which would be too broad to coalesce in tn, become detectable sources of GW with the help 
of a moderate kick velocity. If the companion is already a compact star, the orbital evolution 



is driven exclusively by GW emission (Section A1.4 ). Unless the kick velocity is so big that it 



makes the semimajor axis very large, these binaries will also merge in a time interval shorter 
than the one following from the evolution without a kick. 

These qualitative considerations explain the outcomes of numerical simulations with many 
trial systems. We are interested in results averaged over many systems with different input 
parameters. These results are presented below. As expected, a moderate kick velocity increases, 
on average, the rate of compact star mergings. 

2.3.3 Coalescence rates of compact binaries 

We can now present the results of our numerical calculations for the coalescence rate of compact 



binaries in a typical galaxy |81|. The total mass of a model galaxy is assumed to be IO^^Mq. 
We adopt a constant star formation rate defined by Eq. (|3|). It is believed that Eq. (^) reflects 
well the situation in a galaxy like our own Milky Way. 

In Fig. m we plot the NS-I-NS, BH-I-NS, and BH-I-BH merging rates as functions of the 
kick velocity parameter wq in the distribution (^). The calculations were performed for discrete 
values of wq, but the resulting points are joined by smooth curves. The BH formation parameters 
were taken from the range = 15 — 50Mq with ksH = 0.25. Both the high mass loss and 
the low mass loss stellar winds were considered. The broad range of M^, and the uncertainty 
in the stellar winds have contributed to the spread of the results for BH-I-NS and BH-I-BH 
systems. The NS-I-NS systems arise from relatively low mass stars, so they are less sensitive to 
the uncertainty in the stellar wind. It is seen from Fig. ^ that the NS-I-NS rate lies in the range 
~ 3 X 10~^-~ 3 X 10"^ per year. The rates of BH-I-NS and BH-I-BH mergings are 10-100 times 
lower. For the limiting case of zero kick velocity (wq = 0) our rates agree with the independent 



estimates of Tutukov and Yungelson |82]. In the same limit wq = 0, our rate for NS-I-BH binaries 
(~ 10~^ per year) is smaller than the estimate by Bethe and Brown ||7^, who obtained the rate 
~ 2 X 10~^ per year. However, we believe that their estimate was derived from a somewhat 
simplified picture of binary evolution. 

As expected, the BH-I-NS and BH-I-BH rates have a tendency to grow with the increase 
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Figure 3: NS+NS, BH+NS, and BH+BH merging rates in a 10^^ Mq galaxy as functions of the 
kick velocity parameter wq for Lyne-Lorimer kick velocity distribution (|^. Star formation rate 
in the galaxy is assumed constant. BH formation parameters are M^: = 15-50 Mq, ksH = 0.25. 



of kick velocity from zero. This is seen on the graph as the rise of the NS+BH and BH+BH 
curves for small and moderate wq (up to wq ^ 100 km/s). For much larger values of wq, the 
kick velocity contributes mostly to the disruption of binary systems, and this is why the curves 
have a tendency to turn down. Generally speaking, the NS+NS rate should also grow for small 
deviations of wq from zero. However, since the NS mass is smaller than the BH mass, the 
increase of the NS+NS rate takes place for only a small value of wq, not resolvable on the graph. 
For larger values of wq, the kicks mostly disrupt the binaries, and the NS+NS curve goes down. 
The value of wq preferred by the radio-pulsar observations lies in the range 200-400 km/s. 

For a broad range of used parameters and despite all the remaining uncertainties, the results 
of evolutionary calculations show that the number of coalescing BH+BH pairs is only a factor 
10-100 smaller than the number of coalescing NS+NS pairs. This relationship may have a 
simple explanation and can be traced back to the initial conditions of star formation. The line 
of argument is as follows. Let us take the NS mass at IAMq (a typical mass well confirmed 
by existing observations), and the BH mass at 8.5Mq (the mean value for BH candidates from 
Table |2|). Assume that the lower initial mass of NS progenitors is M{NS) ~ IOMq, while the 
threshold for a BH formation is at the maximum of the estimates quoted above: M{BH) > 
Mcr = SOMq . Applying the Salpeter initial mass function for the formation rate of stars in the 
Galaxy (see Eq. (|)): 

~ lM0(M/M0)-2-35 



dtdM 

and using the lower limits of integration, one finds 

N{M > 8OM0) _ /8OM0 
N{M > IOMq) ~ V IOM0 



-1.35 

~ 0.06. 



This ratio should be valid for binary stars too. It is reasonable to expect then that despite 
differences and complexities of binary evolution, the ratio of coalescence rates will be given, 
approximately, by the same number: 



TZbh /SOMq 



Uns VlOMe 



-1.35 

~ 0.06. (11) 
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This expectation turns out to be in rough agreement with the results of detailed evolutionary 
calculations presented above. 

The derived rates TZq for a single galaxy can be extrapolated to larger volumes. For the 
purposes of GW detection it is important to know the rate of events from distances accessible 
to the instruments in LIGO, VIRO, GEO-600. These are large distances up to and above 
100 Mpc (see Section 0). In such a large volume one can regard galaxies as being distributed 
homogeneously, and at the same time, one can neglect effects of cosmological evolution on star 
formation initial conditions, etc. To derive the average density of Galactic events in a large 
volume one can use different approaches. One possibility is to use the luminosity of galaxies 



per Mpc^ (as in |24|). Alternatively, one can rely on the estimate of density of baryons bound 
in stars. The baryon density is often expressed in terms of the dimensionless parameter 
ilf, = ph/ Per, where Hq is the present value of the Hubble parameter and pcr = 'AHq/S-kG is 
the critical density. Then, one can relate the Galactic rate IZq per a lO^^M© galaxy with the 
volume rate TZy per 1 Mpc'^: 

where e is the fraction of binary stars and /170 = Hq/{70 km/s Mpc). This estimate agrees with 
that of assuming e = 1 (all stars are binaries). The available astronomical measurements 
of the total baryon budget give ~ O.OOlS/iyQ^ in galactic disks and ^If, 0.003/iyQ^ in bulges 
of spirals and ellipticals ||8^ (as well as the somewhat larger values fS^). On the other hand, 
estimates of based on the primordial nucleo-synthesis considerations give as much as Qhh'^ = 
0.016, but this number can also be a factor of 2 smaller [85|. Formula (|l^) can be rewritten as 



TZv = O.ITZg — . (13) 

0.5 0.016 VMpcy ^ ^ 

When comparing our numerical simulations, described below, with qualitative estimates, we rely 
on the relationship 

This result for TZy is based on TZg for spiral galaxies. For elliptical galaxies the star formation 
process is more like an instantaneous event rather than a continuing process described by (^). 
The coalescence rates have been calculated for elliptical galaxies too. However, it was shown [^] 
that the contribution of elliptical galaxies to the coalescence rates from the discussed distances 
is only about 10-20 %. 



3 Detection Rates 

Having found the coalescence rates TZy for binaries of different nature, one can now predict the 
detection rates of these binaries in a given GW detector. We argue that binary black holes have 
a better SNR than of binary neutron stars, and, despite their lower abundance, the BH+BH 
and BH+NS pairs should be seen more often than NS+NS pairs. In the first subsection, we 
derive the detection rates that are based on the TZg described above. In the second subsection, 
we discuss possible modifications to our conclusion in connection with the recently proposed 
scenario |8^, which applies to very massive stars. Since the proposed scenario can affect only 
the BH+BH detection rates, we concentrate on these systems emphasizing the important role 
of kick velocities. 
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3.1 Detection rates in the usual picture 

The rate of NS+NS coalescences is higher than the rate of NS+BH and BH+BH coalescences. 
However, the BH mass is significantly larger than the NS mass. A binary involving one or 
two black holes, placed at the same distance as a NS+NS binary, produces a significantly larger 
amplitude of gravitational waves (see Section | and ^ . With a given sensitivity of the detector 



(fixed SNR), a BH+BH binary can be seen at a greater distance than a NS+NS binary. Hence, 
the registration volume for such bright binaries is significantly larger than the registration volume 
for relatively weak binaries. The detection rate of a given detector depends on the interplay 
between the coalescence rate (spatial density of sources) and the detector's response to sources 
of one or another kind. 

Coalescing binaries emit gravitational wave signals with a well known time-dependence 
(waveform) . This allows one to use the technique of matched filtering Q . The signal-to-noise ra- 
tio S/N depends mostly on the "chirp" mass of the binary system = (Mi+M2)"^/^(MiM2)^/^ 
and its distance r. The accurate formula for S/N is presented in SectioriS (formula (|134| )). Here, 
we will use its simplified version which is sufficient for our purposes (Q, see also fSql ): 

I = 3-^/^vr-^/3^:^/-V6//,_(/) . (15) 

At a fixed level of S/N, the detection volume is proportional to and therefore it is proportional 
to A1^/^ . The detection rate P for binaries of a given class is the product of their coalescence 
rate TZy with the detector's registration volume oc A4^/^ for these binaries. 
Let us start from a qualitative discussion of the expected ratio 

(16) 



T^NS ^NS \Mns 

where Vbh and Vj^s refer to BH+BH and NS+NS pairs, respectively. Here, we discuss the 
ratio of the detection rates, rather than their absolute values. The derivation of absolute values 
require detailed evolutionary calculations which will be discussed later. As a rough estimate for 



T^BH /T^NS one can take Eq. (pTl). Then, Eq. ([Lq) gives a remarkable result: 



Vns \WMqJ \1A0Mq 

This ratio becomes even larger than 5.5, if one takes Mcr < 80Mq as usually assumed. Thus, 
the registration rate of BH mergers is expected to be higher than that of NS mergers. This 
estimate is, of course, very rough, but it can serve as an indication of what one can expect from 
detailed calculations. 

In Fig. ^ we display the results of numerical calculations for the absolute registration rates 
of various binaries. The detector sensitivity is taken as hrms = W~^^ at / = WOHz, as expected 
for initial instruments in LIGO, VIRGO, GEO600. It is assumed that S/N = 1. Since the most 
interesting results refer to systems with black holes, we vary the black hole formation parameter 
^BH- The calculations were performed assuming the Lyne-Lorimer kick velocity distribution 
with wq = 400 km/s. The vertical dispersion of the results is due to the uncertainty in the 
parameter Mer- 
it is seen from the graph that under the formulated conditions one can expect to see a couple 
of NS+NS coalescences in 1-3 years of observations, at a S/N = 1. These systems are located, 
roughly, at a distance of 100 Mpc. The SNR is higher for closer systems, but the expected 
event rate would be lower; for more distant systems, the event rate increases, but the SNR 
becomes smaller than 1. So, it is unlikely that NS+NS coalescences will be detected by the 
initial instruments. 
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Figure 4: The detection rate V of GW events in a detector with the sensitivity hrms = W at 
frequency 100 Hz and the signal-to- noise level S/N = 1, as a function of BH formation parameter 
ksH- The calculations were performed for the Lync-Lorimer kick velocity distribution with 
wq = 400 km/s. The spread of T> at fixed ksH is due to variation of the parameter Mcr from 
15 Mq to 50 Mq. The bottom rectangular area is drawn for binary NS coalescences. Their 
rate is independent of ksH (as it should be) and predicts a couple of events per 1-3 years at 
this level. The total detection rate can be 2-3 orders of magnitude higher then the NS+NS 
rate at the expense of BH+BH and BH+NS coalescences. The hatched area shows the region 
of the most probable parameters for the low stellar wind mass loss scenario. Inside this region 
the outcomes of calculations are in agreement with the upper limit on the galactic number of 
binary BH with radiopulsars (less than 1 per 700 single pulsars) and the galactic number of BH 
candidates similar to Cyg X-1 (from 1 to 10). 
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The situation is significantly better for systems involving black holes. As is seen from Fig. ^, 
the total registration rate of all binaries, including BH+BH and BH+NS pairs, can be 2-3 orders 
of magnitude higher than the registration rate of NS+NS systems alone, mostly at the expense 
of massive BH pairs. This is true unless the ksH parameter is very small, ksH < 0.4. The 
hatched area shows the results of calculations with the stellar wind parameters taken from the 
"most probable" region. This means that under this choice of the parameters, the outcomes of 
other evolutionary tracks are in agreement with observations, namely, in agreement with the 
upper limit on the number of binary pulsars with BH (less than 1 per 700 single radio-pulsars) 
and with the number of Cyg X-l-like BH candidates (from 1 to 10 per Galaxy). Inside this 
region, one should count on 100 registrations (at the level S/N = 1), mostly from BH mergers. 
The mean total mass of the BH pairs in the hatched area is around M = SOMq. The simplified 
formula (^), used in the construction of Fig. ^, overestimates the S/N for pairs heavier than 
3OM0, as shown in Section ^. However, a correction for more massive binaries is not expected 
to change significantly the derived total registration rate. 

For a reliable detection, the S/N ratio should be at least 2 in each of a network of four or 
more antennas. Then, the calculated detection rates should be decreased by at least a factor 
{S/N)^ = 8. This is because of the scaling S/N oc 1/r and P oc cx {S/N)~^. Then, 
the expected detection rate of merging BH+BH pairs is up to 10 events per year. As will be 
explained in Section^, the S/N ratio is somewhat different for the three different instruments: 
LIGO, VIRGO and GEO. For a coalescing pair with a total mass M = SOM© the SNR is roughly 
4 for sources at a distance of 100 Mpc. (The performance of VIRGO is expected to be better 
than that of other instruments, since VIRGO will be more sensitive at lower frequencies and can 
track the binary for a larger number of cycles.) If one is satisfied with S/N = 2, the accessible 
radius increases to r = 200 Mpc. Then, the calculated detection rate (several per year) is in 



agreement with formula (14) if one takes for coalescing black holes a reasonable galactic rate 
TZg = 3 X 10^^ and r = 200 Mpc. In its turn, this value for TZq fits well the event rate derived 
from numerical simulations, as displayed in Fig. |3[ 

Thus, taking into account all the remaining uncertainties, we conclude that the initial net- 
work is likely to see each year 2-3 coalescing black hole binaries with the total mass around 
3OM0, at an SNR level of about 2-3. 

3.2 Non-standard scenarios and effects of kick velocities 
on BH+BH detection rate 



Some of the recent evolutionary calculations [S7| assume that the primary stars with initial 



masses Mi > 4OM0 never fill their Roche lobes, so that the components of the binaries evolve 
like single stars. As a result, the binary BH systems would be too wide to merge in tn- Although 
we think the scenario |87| will face observational difficulties, since it will lead to the too small 



number of binaries involving a BH and a massive blue star (Cyg X-l-like systems), we consider 
it necessary to follow in detail the possible fate of binary BH systems. We argue that the 
kick velocity accompanying the BH formation increases the eccentricity of the binary, decreases 
its coalescence time, and thus keeps the detection rate at almost the same level as discussed 
in section p.l| . In addition, the kick velocity leads to interesting modifications in the relative 
orientations of the black hole spins with respect to each other and with respect to the orbital 
angular momentum. 



We have adopted the proposed scenario and have carried out population synthesis cal- 



culations by varying the kick velocity parameter. The binary BH merging rate was derived for 
a model galaxy of 10^^ Mq (assuming that all stars are formed in binaries) with a constant star 
formation rate. For simplicity, the kick velocity distribution was taken as a delta-function. The 
more complicated distributions do not change the results significantly and are not commented 
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Figure 5: New scenario — superhigh wind. Left: BH+BH merging rate calculated for a 10 Mq 
galaxy with a constant star formation rate, as a function of the kick velocity during BH formation 
with M(.r = 35Mq, for k^H = 0.5 and 0.75. Right: Detection rate of BH+BH mergings by the 
initial laser interferometers (/irms = 10~^^ at / = 100 Hz), as a function of the kick velocity 
during BH formation. 

upon here. The results of these calculations are shown in Fig. ^. The left panel shows the 
merging rate, while the right panel shows the detection rate. The detection rate of binary BH 
coalescences is given for initial laser interferometers (/irms = 10~^^ at / = 100 Hz) as a function 
of BH kick velocity. It is seen from the plot that the merging rate and the detection rates in- 
crease rapidly with the kick. The merging and detection rates reach the maxima TZ ~ 2.5 x 10~^ 
yr~^ and D ~ 20 detections per year for w ~ 120 km/s. Since T> ~ M.^^'^TZ, the lZ{w) and Viw) 
functions have similar shapes. 

Obviously, the kick velocity imparted to newly born black holes makes the orbits of survived 
systems highly eccentric. It is important to stress that some fraction of binary BH retain their 
large eccentricities up to the late stages of their coalescence. This signature should be reflected 
in their emitted waveforms and should be modeled in templates. 

The asymmetric explosions accompanied by a kick velocity change the space orientation of 
the orbital angular momentum. On the other hand, the star's spin axis remains fixed (unless the 
kick was non-central). As a result, some distribution of the angle between the BH spins and the 
orbital angular momentum (denoted by J) will be established [^]. It is interesting that even 
for small kicks of a few tens of km/s an appreciable fraction (30-50%) of the merging binary BH 
should have cos J < 0. This means that in these binaries the orbital angular momentum vector 
is oriented almost oppositely to the black hole spins. This is one more signature of imparted 
kicks that can be tested observationally. These effects are also discussed in a recent paper ||90t| . 

Thus, to conclude this analysis, we stress again that binary black hole coalescences remain 
the most likely sources to be detected first by the initial network of laser interferometers. 

4 Transients and Continuous Gravitational Waves 

In this Section we will discuss two distinct types of signals: (1) transient events, that last few to 
several milliseconds, which, on astronomical grounds, are expected to occur but emit waves of 
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unknown phase evolution as in the case of supernovae and (2) continuous radiation, that last for 
several days or longer, from either newly born neutron stars or old recycled neutron stars. The 
strengths, duration and shapes of these signals is rather speculative and highly uncertain. On 
general physical grounds we should expect such sources to exist and every effort should be made 
in searching for these sources by taking the best advantage of current knowledge. However, 
since the astrophysical uncertainties are so large, we shall keep the discussion of this topic quite 
qualitative. 



4.1 Transients 

4.1.1 Supernovae and asymmetric explosions 

Supernovae (of type II) are associated with violent mass ejection with velocities of order 0.001c 
and the formation of a compact remnant — a neutron star or black hole. The event has at its 
disposal the difference in the gravitational binding energy of the pre-collapse star and the newly 
formed compact star which, neglecting the former, is: 

\E\^3x 10^3 ( -JL_] erg. (18) 

99% of this energy is carried away by neutrinos, about 1% is transferred as kinetic energy of 
ejecta, a fraction 10"'* of the total energy is emitted in the form of electromagnetic radiation. 
Depending on how asymmetric the collapse is, some fraction of the total energy should be 
deposited into gravitational waves; spherically symmetric collapse, of course, cannot emit any 



radiation. According to numerical simulations (see |91] for a review) one might expect up to 
10-"^ of the total energy to be emitted in gravitational waves. Together with uncertainties in the 
event rate, this is not a very encouraging prognosis for the initial instruments ^. It appears 
that a star collapsing to form a black hole is also not particularly well suited for detection by the 
existing resonant detectors and forthcoming interferometers [^2|. However, second generation 
interferometers should be able to see a supernova event as far as the Virgo super-cluster which 
contains about 200 bright galaxies and at least twice as many faint galaxies. In addition, there 
are a few other smaller clusters within that distance as also a large number of field galaxies. 
Therefore, the supernova event rate for these instruments could be as large as tens per year. 
Such observations would undoubtedly be of great interest and would shed light on hitherto 
un-understood processes that occur when a star collapses to form a compact object. 

An asymmetric collapse generates gravitational waves but is not necessarily accompanied 
by a change of linear momentum of the exploding star. However, as we emphasized above, the 
observations of single and binary radio-pulsars require the presence of a significant kick velocity 
imparted to a newly born neutron star. The kick velocity can be as large as 200-400 km/s. This 
testifies for a non-axisymmetric explosion leading to an additional linear momentum acquired by 
the neutron star. A possible reason for this asymmetry can be related with the asymmetry of the 
neutrino emission during the collapse [^3|. Using the observed kick velocities one can evaluate 
asymmetry of the explosion and calculate the amplitude of emitted gravitational waves. Nazin 



and Postnov [94| estimate that the mean energy carried away by the gravitational wave burst 
can reach the value Eqw = 5 x 1O~^M0C^. It is interesting to note that the radiated signal 
belongs to the category of bursts "with memory" [ p5| , |96[ . However, the estimated strength of 
the asymmetric explosions is still too low to count on them as reliable sources for first detection. 



4.1.2 Bar-mode and convective instabilities 

Some amount of gravitational waves may be emitted during two special stages of the collapse: 
(1) rotation-induced bars and (2) convective instabilities set up in the core of the newly born 
neutron star. 
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If the core's rotation is high enough it may cause the core to flatten before it reaches nuclear 
density, leading to an instability that transforms the flattened core into a bar-like configuration 
which spins about its transverse axis. Some of these instabilities could also fragment the core 
into two or more pieces which then rotate about each other. Both are efficient ways of losing 
energy in the form of gravitational waves. It is estimated [97| that the waves could carry up to 
IO'^Mqc^ in a few ms. LIGO and VIRGO detectors can see such an event at a distance up to 
50 Mpc, or about 5-10 Mpc if the waves come off at 1 kHz. GEO will also be able to see these 
events provided the signal comes off around 200 Hz |98|. 

Convective instabilities in the core of a newly born neutron star, which last for about a 
second after the collapse, are likely to produce gravitational radiation due to anisotropic mass 
distribution and motion . Miiller and Janka |100| ] find that the gravitational wave amplitude 
of a source at 100 kpc is /i ~ 10~^^ and that the waves would come off at about 100 Hz. Since 
there would be about 100 cycles, one can enhance the amplitude to about 10"^'^ provided we 
know the development of the signal but this is still far too weak to be detected beyond about 
10 kpc with high confidence. 



4.1.3 Merger vi^aves 

The physics of the merger phase that follows the adiabatic inspiral regime of a compact binary 
coalescence is not known. Presently, it is a subject of active research and contains a lot of 
uncertainties. Some authors expect that this phase could result in detectable amounts of ra- 
diation [^] while others conclude that the merger phase is just a continuation of the inspiral 
phase adding on the order of 1 more cycle to the inspiral [101]. There is a lot of effort to solve 



this important problem by using semi- analytical and numerical techniques. It is likely that a 
solution will be in place by the time the first detectors will begin to operate. The fact that the 
coalescence waves will be preceded by inspiral waves makes the search easier, though it is not 
inconceivable that while the former may be observable with the aid of accurate search templates, 
the latter may not be. Flanagan & Hughes |^8[ estimate that for binary systems of total mass in 
excess of 25M0, coalescence waves are likely to be significantly stronger than the inspiral waves. 
The span of a detector is larger for heavier binaries and therefore, these authors conclude, it is 
likely that the first gravitational wave events will be the merger phases of massive binaries. 



4.1.4 Sub-stellar mass black hole binaries 

The results of recent micro-lensing experiments J102|] ha ve revealed massive compact halo objects 
(MACHOs) of mass O.btHMQ. Nakamura et al., |103|] argue that if MACHOs are black holes 
then they must have formed in the early Universe and they estimate that our Galaxy may 
contain about 10^ black hole binaries with inspiral time scales less than the Hubble time. If this 
is the case then the rate of MACHO coalescences in our Galaxy is ~ 5 x 10~^ yr~^, implying an 
event rate of few per year within 15 Mpc. As we shall see in the Data Analysis Sections, first 
generation interferometers should be able to detect the final inspiral phase of these systems. 



4.1.5 Quasi-normal modes 

A compact binary coalescing as a result of gravitational radiation reaction would most likely 
result in the formation of a single black hole. The newly formed hole will initially be some- 
what non-spherical, and this dynamical non-sphericity will be radiated away in the form of 
gravitational waves. The late time behaviour of this radiation is well studied in the black hole 
literature and there are detailed calculations of the (quasi) normal modes for both static, i.e. 
Schwarzschild, and stationary, i.e. Kerr, black holes. In all cases the time-evolution of the 
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emitted radiation is well-modeled by a quasi-periodic signal of the form 



h{t; T, uj) = Ae'^l'^ cos{ujt) (19) 

where r is the decay time-scale of the mode in question and lo is the angular frequency of the 
mode, both of which depend on the black hole mass and angular momentum. In all but the 
extreme Kerr black holes (extreme Kerr black holes are those that are spinning at the maximum 
possible rate) the only dominant mode, i.e. the mode for which the decay time is the longest 
and the amplitude is the highest, is the fundamental mode whose frequency is related to the 
mass M and spin a of the black hole via uj = [1 - 0.63(1 - a)0-3]/M where M is the mass of 
the black hole in units G = c = 1 and a = J/M is the spin angular momentum of the hole in 



units of black hole mass. The decay time r is given by r = 4/[lij(1 — a)'^''^^]. (See Ref. [104] and 
references therein for details.) 



It is estimated [38[ that during the quasi-normal mode ringing of a black hole the energy 



emitted might be as large as 3% of the system's total mass. By matched filtering (cf. Sec. 8.1 
it should be possible to detect quasi-normal modes, in initial interferometers, from black holes 
of mass in the range GO-IO^Mq and at a distance of 200 Mpc. Binary black hole mergers should 
result in the emission of such a ring down signal during the late stages. Thus, inspiral signals 
emitted before the merger might aid in identifying the quasi-normal modes. 



4.2 Continuous waves 

Our Galaxy is expected to have at least 10^ spinning neutron stars that form roughly at a rate of 
one every 30 years. Some population of neutron stars is in binaries. There are a number of ways 
in which a single spinning neutron star could radiate away gravitational waves (if the neutron 
star is axisymmetric, of course, then there will be no gravitational wave emission): (1) Neutron 
stars normally spin at high rates (several to 500 Hz) and this must induce some equatorial 
bulge and flattening of the poles. The presence of a magnetic field may cause the star to spin 
about an axis that is different from the symmetry axis leading to a time-varying quadrupole 
moment. (2) The star may have some density inhomogeneities in the core/crust set up during 
its formation and/or subsequent convectively unstable motions of the core. (3) The presence of 
an accretion disc, with its angular momentum not necessarily aligned with that of the neutron 
star, can potentially alter axisymmetry. That and electromagnetic radiation reaction torques 
can induce and sustain wobble. (4) The normal modes of the neutron star fluid (radial and 
other oscillations) can extract rotational energy and re-emit in the form of gravitational waves. 
(5) There are certain classical and relativistic instabilities in the neutron star fluid which may 
cause the star to radiate away energy in the form of gravitational radiation. In what follows we 
will only discuss a sample of recent work on the radiation from spinning neutron stars. 



GW amplitude from spinning asymmetric neutron stars If Izz is the moment of inertia 
about the spin axis of a neutron star emitting gravitational waves at a frequency / then the 
gravitational amplitude at a distance r is: 

where e is the ellipticity of the star. In a simple model of an equatorial plane of elliptical cross 
section of semi-major axis oi and semi-minor axis 02, the ellipticity is e = 1 — 02/01- The 
ellipticity is an unknown but one can obtain an upper limit on it by attributing the observed 
spin-down of pulsars P to gravitational radiation back reaction, namely that the change in the 
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rotational energy E = 70^/2 is equal to gravitational wave luminosity. Then, the ellipticity is 
related to the spin-down rate of a pulsar via 

P \ 3/2 / p \ 
5.7x10-Mt.-^ — ^ . (21) 



10-2 sy v^^~^^/ 

Since one knows the observed values of P and P one can obtain an upper limit on e using the 
above equation. Following this method one can find that for the Crab pulsar e < 7 x 10"^ and 
the gravitational amplitude is h < 10"^^. One concludes that Crab will have a sufficiently large 
amplitude to be observable with the aid of GEO detector in 1 year of continuous observing, 
if all of its spin down can be attributed to the emission of GW. It is unlikely that the ellipticity 
is so large. Yet, the prospect of seeing Crab event at a 10th or a hundredth of this ellipticity is 
quite good with first /second generation of interferometers. 

In the next two paragraphs we will discuss some new developments in relativistic astrophysics 
that could lead to potential GW sources. 



Relativistic instabilities in young neutrons stars Chandrasekhar |105] and Friedman & 



Schutz |10f:] discovered an instability (now called the CFS instability) in the fundamental, or "f", 
mode of a neutron star fluid, arising as a result of GW emission. The mode goes unstable above 
a critical spin frequency of the star and progressively grows, instead of decaying, by emitting 
gravitational waves. 

The physics behind this instability can be understood in the following manner: Imagine 
exciting a mass-quadrupole mode - that is a non-uniform distribution of mass - in a non- 
spinning star. The mass inhomogeneity will travel on the surface of the star and this dynamical 
asymmetry will cause the star to radiate gravitational waves. After a while, all the energy in 
the mode will be radiated away and the mode will decay. Now consider a spinning neutron 
star in which a co-rotating and counter-rotating modes are excited. These fluid modes have a 
certain pattern speed on the surface of the star. For low spin rates both these modes will decay 
in course of time by emitting gravitational waves. But as the neutron star is spun up above a 
critical rate, to an external inertial observer both modes will appear to be traveling in the same 
sense as the rotation of the star. Therefore, the mode counter-rotating relative to the star will 
emit positive angular momentum, causing the angular momentum associated with the mode to 
enhance, or for the amplitude of the mode to increase. In other words, a mode counter-rotating 
relative to the star, but seen co-rotating relative to the inertial observer, can only emit negative 
angular momentum which causes its own angular momentum to increase. The energy for this 
enhancement is supplied by the spin angular momentum of the neutron star. Thus, while the 
mode co-rotating with the star's spin will decay, the mode counter-rotating with the spin will 
grow in amplitude and emit more and more radiation. This will go on until the mode has sucked 
out enough angular momentum of the star to make the counter-rotating mode appear to be 
counter-rotating with respect to an inertial observer too. 

It is suspected that the CFS instability will not work in the presence of viscosity and hence 
it may be unimportant in old neutron stars. However, newly born neutron stars will be very 



hot and viscous forces may be insignificant in them. Recently, Andersson [107] discovered 
another class of modes called r-modes, which — unlike CFS modes that are mass-quadrupole 
moments — are current-quadrupole moments, that are unstable at all spin frequencies. The 
role of these modes may be significant in young neutron stars. It is proposed that r-modes are 
responsible for the limit on the spin frequencies of newly born neutron stars [|108(| . Owen et al. 
|109| ] have computed the efficiency with which these modes extract energy out of the system 
and the expected gravitational wave amplitude from isolated neutron stars, as well as from the 
ensemble of all sources up to cosmological distances. They conclude that the second generation of 
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interferometric antennas will be able to distinguish such a background by coincident observations 
with a nearby bar detector. Such observations should provide independent knowledge of the 
distribution of galaxies in the high-red shift Universe as also on the star formation history via 
the observation of r-modes associated with the formation of neutron stars and pulsars. 

Neutron stars in X-ray binaries In the recent years, Rossi satellite observations of the 
X-ray emitting binaries have shown high-frequency quasi-periodic oscillations (QPO) in their 
X-ray power spectra. Some authors believe that these QPO could be a result of beating of 
two frequencies one of which is that of the neutron star. The neutron star spins inferred in 
this way seem to lie in a narrow range of 250-350 Hz and are all within 20 % of 300 Hz. A 
neutron star may be born with a high spin rate (several 100 Hz) but quickly spins down to 
moderate rates (several 10 Hz). In a binary system, when the companion becomes a red giant 
the neutron star starts accreting mass and angular momentum. Though the mass accretion 
rate is very low (M ~ 10"^*^ Mq yr~^) the accretion of angular momentum can spin up a 
neutron star and heat its crust substantially. Indeed, the millisecond pulsars are believed to 
be old pulsars in binaries recycled in this way. It is puzzling as to why the spin frequencies of 
neutron stars in X-ray binaries are all in a narrow range. Bildsten suggests |110| that absence of 



efficient heat transport processes make it possible to set up temperature gradients in accreting 
neutron stars. Provided that the temperature distribution has large scale asymmetry then 
the temperature sensitive electron captures in the deep crust can build up the mass quadrupole 
(~ 10~ ''Mi?^) needed to radiate away accreted angular momentum and limit the spin frequency. 
This mechanism is present only during accretion and decreases rapidly once the accretion halts. 
The frequency of GW radiation will be known in advance since one knows the spin frequency of 
the pulsar via X-ray observations. It is argued |110| ] that the gravitational wave strength will be 
h ~ (0.5-3) X 10~^^ for many of these sources and that the LIGO/VIRGO and signal recycled 
GEO600, can detect the strongest of these sources, Sco X-1, at an SNR of 5 with a few years of 
integration. 

5 Astrophysical stochastic backgrounds of gravitational waves 

Coalescing binaries of compact stars are at the center of our attention. A gravitational wave 
signal will be monitored by a detector as long as the changing frequency of the source sweeps 
through the detector's window of sensitivity. The detection of a useful gravitational wave source 
is limited by the instrumental noise and by a possible gravitational wave noise produced by 
other sources. If a large population of astrophysical sources emit overlapping gravitational 
waves, the resulting signal will be perceived by the detector as a gravitational wave noise. In 
principle, this noise could be larger than the instrumental noise. The astrophysical backgrounds 
of gravitational waves are important in their own right, but we are discussing them mostly 
as possible noises preventing the detection of a signal of a greater scientific importance. The 
aim of our discussion is the derivation of the gravitational mean square amplitude in a given 
frequency interval A/ and comparison of this number with the instrumental noise in the same 
interval. If the background is below the instrumental noise, it will not prevent the detection of 



the useful signal. In subsection 5.1 we derive general formulas for the performance of a large 
number of unresolved sources. In subsection ^.1.1| we discuss the population of rotating neutron 
stars in our Galaxy. They could constitute a major danger for the ground-based interferometers. 
The analysis shows that the danger arises only under quite unrealistic assumptions about the 



parameters of this population. In subsection 5.1.2 we discuss the stochastic gravitational wave 



background produced by binary white dwarfs in the Galaxy. This background dominates the 
LISA instrumental noise from ~ 10~^ Hz up to ~ 10~^ Hz, but leaves the gravitational wave sky 
transparent at lower and higher frequencies. Finally, subsection |5.2| comments on astrophysical 
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backgrounds of extra-Galactic origin. Generically, these backgrounds are one order of magnitude 
lower than those produced by the Galactic sources. 



5.1 Unresolved sources in our Galaxy 

A large collection of independent sources produces signals whose intensities add. Consider 
identical sources located at approximately the same distance r from the observer. The resulting 
gravitational wave field is characterized by the r.m.s. amplitude h^'- 



hN = hiVN, (22) 

where hi is the averaged amplitude of a single source. We are interested in the narrow frequency 
interval from / to / + A/. The radiating sources gradually change their frequency and pass 
through the window of interest. The crossing of the window can occur either on the way from 
lower to higher frequencies, as in the case of binaries, or in the opposite direction, as in the case 
of rotating neutron stars. To find we need to know the number N{f, A/) of the radiating 
systems in the discussed frequency interval. Denote by TZ the rate at which the sources appear 
in the window. A source is present in the window during the time At, where At = A///. Thus, 

A(/,A/)=7^At = 7^^, (23) 

and 

I aT 

(24) 




If the frequency evolution is driven by gravitational waves only, the quantity / is determined 
by the gravitational radiation damping. For example, in the case of a binary system, the / is 
given by Eq. ( AlSj ). Alternatively, the quantity / can be determined by the electromagnetic 



radiation damping. This takes place in the case of highly magnetized rotating neutron stars 
whose / is determined by electromagnetic, rather than gravitational, losses. The appearance 
rate TZ of sources of a given population in the discussed window Af is equal to the birthrate or 
to the coalescence rate of the sources as such (see, for example, T^coaies in Sec. ^). This is true 
as long as the number of sources remains constant during the frequency evolution throughout 
the window, what we always assume. 

For sources driven by gravitational radiation only, formula (^) can be expressed in terms 
of the energy of a single radiating system and the rate TZ. One needs to use Eq. ( A16| ) in order 



to express hf in terms of dE/dt, and to take into account the relationship dE/dt = {dE/df)f. 
Then, Eq. ( p4|) takes the form 

,2 G n dE^^ 



Usually, the energy of the radiating system is a power-law function of the frequency /: E{f) 
For example, in the case of a binary star in circular orbit one derives from Eq. (A6): 



E{f) = ^M''\^ffl\ (26) 

where M is the chirp mass. So, in the case of binaries, a = 2/3. For a non-axisymmetric 
rotating star E{f) ~ Z'^, so that a = 2. Thus, Eq. (p5|) takes the universal form 



2 _ G IZ „^,,,A/ 



= -^r^^Eif)-^. (27) 
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The quantity hj^/Af is the mean square noise amphtude h'j (with dimensionahty Hz ^] 
which participates in the expression 



J max 

{h')= J h}df, (28) 



/mill 

and can now be compared with the frequency-dependent instrumental noise. We will work with 
the dimensionless spectral amplitude 




= (29) 

Obviously, the independent sources in a frequency bin A/ form a kind of stochastic back- 
ground if the number N(f,Af) is much larger than 1. A source becomes resolvable, if this 
number is of order 1. For a collection of evolving sources one can find the limiting frequency /um 
at which this happens. For fixed A/ and TZ, one uses the concrete function /(/) (arising due to 
the gravitational reaction force or by some other reasons) and finds /um from the requirement 

= 1. (30) 

/ 

For a collection of binary stars one uses Eq. ( AlSp and finds 



/h. -(1.2 X 10-3Hz)7^^/- f ^..r^T' ( J^)"'" • (31) 



3xlO-8Hzy V0.52M, 



The appearance (coalescence) rate TZ^qq is chosen for compact white dwarfs, which are expected 
to coalesce in our Galaxy once per 300 years. The chirp mass is normalized to 0.52 Mq which is 
true for two CO white dwarfs with masses 0.6 Mq. This estimate will be needed in subsection 



5.1.2 which discusses the LISA noise. 



5.1.1 Noise from old neutron stars at frequencies of ground-based interferometers 

Rotating neutron stars as sources of gravitational radiation can be roughly divided into two 
populations: One consists of old neutron stars with relatively weak magnetic fields and small 
electromagnetic losses. Their rotational frequency slowly decreases due to gravitational wave 
damping. Another population consists of young highly magnetized neutron stars. Their ro- 
tational frequency decreases much faster due to electromagnetic damping. In course of their 
frequency evolution, members of both populations cross the window of sensitivity of ground- 
based interferometers, descending from 10^ Hz to 10 Hz. The number of sources simultaneously 
radiating in a given frequency interval is proportional to the birthrate TZ of the population and 
inversely proportional to the velocity of the population flow / through the window. We start 
from old neutron stars and then discuss young neutron stars. 

The frequency evolution of old, rotating, deformed neutron stars is governed by gravitational 
radiation damping. Formula for /, analogous to Eq. ( |A1S| ) for double stars, is 

/ = '-^I^'f^ (32) 

where I is the relevant moment of inertia and e is the ellipticity (deformation) parameter. The 
birthrate TZ of neutron stars in this population can be estimated using the observed fraction of 
millisecond pulsars (which are thought to be old neutron stars with low magnetic fields that have 
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been spun up by accretion in a binary system) among radio pulsars: Nms/^PSR ^ 20/2000 = 
1/100. Remembering that the hfe-time of a mihisecond pulsar is tms = 10^ yr and that the 
life-time of an ordinary radio pulsar is tpsR = 10^ yr, and adopting the Galactic birth rate of 
radio pulsars 1 per 30 years, as for the core collapse supernovae, one can estimate the birth rate 
of millisecond pulsars: 

^ ^ ^ms tpSR ^ „ in-6 -1 

T^ms = T^PSR-Tj : PS 3 X 10 yr 



This estimate is in agreement with the one derived in Ref. [111| from the observed space density 
of millisecond pulsars. 



From Eq. (30) one finds the limiting frequency 



/lim « 53 Hz ( —^—-V^' ( .^,tL.r ) (rrT^J—-^) i-^) ■ (33) 



3 X 10^ yr y V 3 x lO-^Hz / V lO^^g cm2 / V lO^^ 

Thus, at frequencies below 53 Hz the population of old neutron stars is likely to produce a 
stochastic background. Taking E{f) = 7r^//^/2 one derives from Eq. ([29|): 

,«(/) 4,/!^ . 2 . 10-- (^) f ^) . (34) 



V r /V3xl0^yr/ ylO^^gcm^ 

It is interesting that this quantity does not depend on the deformation parameter e, as soon as 
e 7^ 0. The /if and / are both proportional to e^, so that cancels out in the expression for h?^. 
The quantity hN{f) is also independent of frequency /. The numerical level of /iAr(/) is much 
lower than the instrumental noise of initial and advanced ground-based interferometers. The 
/i7v(/) can be increased by two orders of magnitude, and hence the gravitational wave noise from 
old neutron stars becomes marginally detectable, only under the condition that one postulates 
a significantly larger (and, we believe, unrealistic) birthrate TZ for old neutron stars (compare 

with dm 



At frequencies higher than the limiting frequency (^3|), the sources are resolvable during a 
1 year interval of observations. However, to monitor a single neutron star one needs to know 
its exact location on the sky and to take care of the Doppler frequency modulation due to the 
Earth's motion around Sun. 

The young neutron stars differ from old neutron stars in that their electromagnetic energy 

loss 

E,^ = ^^'f\ (35) 
where is the NS magnetic moment, is significantly larger than the gravitational wave loss 

4- = ^^I'e'f. (36) 



The ratio x = Eem/Egw is 



where /Uso = /[//(lO^'^G cm^), e_6 = e/10^^ and 145 = I/lO^^g cm^. For typical parameters 
of young neutron stars one has x ^ 1. The ratio x becomes comparable with 1 only for 
relatively weak magnetic fields, such that the magnetic moment /i satisfies the condition ^ < 
1.5 X 10^^(G cm3)e_6(//100 Hz). The frequency change / is determined by the electromagnetic 
loss and reads 

/=3^^- (38) 
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For fiim one derives 



= 0.5 Hz (7^3o)^/^ ( o ti,^ V^' (/U3o)~'/'(/45)'/=^. (39) 



^3 X 10-8Hz 

The averaged amplitude hi of a single neutron star amounts to 



(40) 

Then, the h]\f{f) found from Eq. ( p^ ) with the help of Eq. ( p8[ ) becomes 

^^(/) = - 3 X 10-^ (^) nli'ire^,,,,^f/100Hz. (41) 

Thus, the gravitational wave noise from young neutron stars is at the same numerical level as 
the noise from old neutron stars. The appearance rate of young neutron stars is much higher 
than that of old neutron stars, but they flow through the window of sensitivity much faster, so 
that there aren't a sufficiently large number of sources to produce a background of a high enough 
level. We conclude, that the populations of neither the old nor new neutron stars do present 
any danger for the sensitivity curves of initial and advanced ground-based instruments. It is 
important that the calculated noise level is below the expected level of the relic gravitational 
waves (see Section 6). 

5.1.2 Noise from Galactic binary white dwarfs in LISA 

The LISA frequency range, 10~^~10~^ Hz, can be contaminated by the gravitational wave noise 
from coalescing binary white dwarfs (WD) and binary neutron stars (NS). The appearance rates 
TZ of the two populations in the sensitivity window are numerically equal to their coalescence 
rates. The binary WD systems are much more numerous than the binary NS systems, and the 
coalescence rate of the former population is significantly higher than that of the latter. The 
binary WD coalescence rate is about 1 per 300 yrs, while the coalescence rate of binary NS is 
about 1 per several 10000 years. At the same time, the chirp masses, Ai ~ O.52M0 for a WD 
binary and A4 ~ 1.22Mq for a NS binary, are not so significantly different. This is why the WD 
background is more important than the NS background, and we consider only the former. 

For a collection of sources consisting of binary stars, one can use Eqs. ( |A13|) , ([aT9|) in Eq. 



(24), or, alternatively, Eq. ( |26| ) and a = 2/3 in Eq. (27). By either way one obtains 

(42) 



10" 



2oAOkpc\/ n y/V / \-2/3/ M 



V r j UoOyrj UO-^Hzj VO.52M0 



Roughly, this is a result of performance of ~ 10^ binaries in the frequency bin A/ = / = 10 
Hz with the averaged amplitude (see Eq. ([A13 )) 



It is necessary to note that Eq.(42) gives the estimate for the amplitude averaged over the whole 



sky, while the real background is strongly concentrated toward the Galactic plane [36|. The 
response of a space-based interferometer should be modulated while the instrument is turning 
in its orbit. In principle, this distinctive feature of the Galactic background can be used in order 



to distinguish it from the backgrounds of cosmological origin [113] 
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In formula (42), the quantity TZ replaces all the astrophysical uncertainties in the binary 
WD evolution. At frequencies higher than 3 x 10~^ Hz the evolution of the vast majority of 
binary white dwarfs is totally controlled by GW emission. So, at frequencies of interest, the 
GW noise is fully determined by the Galactic rate of binary WD mergers and is independent of 
complicated details of binary evolution at lower frequencies. (For examples of calculated spectra 



at all frequencies see [114, |3^, ^4|, |115[| .) The coalescence rate of close binary WD is known only 



up to a factor of few. One way to estimate TZ is based on the search for nearby WD binaries. A 



recent study [116] revealed a larger number of such systems than had been previously believed 
to exist. However, the statistics of such binaries in the Galaxy remains very poor. If coalescing 
binary WD are associated with SN la explosions, as proposed by [|ll7f[ and further investigated 



by many authors (for a recent review of SN la progenitors see [118|), their coalescence rate can 



be constrained using the much more representative SN la statistics. The authors in Ref. [ llg [ 
have concluded that the coalescing CO-CO binary WD remain the most plausible candidates 
responsible for the SN la events. The Galactic rate of SN la is estimated to be 4 x 10~^ per 
year |119| , |120| , which is close to the calculated rate of CO-CO coalescences ~ (1-3) x 10~^. 
The coalescence rate of the He-CO and He-He WD pairs (other possible progenitors of SN la) 



falls ten times short of that for CO-CO WD [118[. As SN la explosions may well be triggered 



by other mechanisms as well, we conclude that the observed SN la rate provides an upper limit 
to the double WD merger rate, regardless of the evolutionary considerations. 



In Fig. 13 we plot the LISA sensitivity curve (thick solid line) calculated for the frequency 
bins A/ = 3 X 10"^ Hz as a function of frequency. The binary confusion limit is shown with the 
dash-dotted line. At frequencies below ~ 4 x 10~^ Hz the binary GW background is produced not 
only by coalescing WD, but by other binaries as well. So, in this part of the graph we rely upon 



numerical calculations [^, 121]. At the limiting frequency ~ 10~^ Hz, individual Galactic WD 
binaries become resolvable in a 1 year observation time, and the binary WD noise drops below 
the LISA sensitivity. It then continues as a noise produced by isotropic distribution of extra- 
Galactic binaries (see Sec. 5.2). Fig. also shows the expected background of relic gravitational 
waves (see Section |6|) . Keeping in mind that the real noise caused by merging Galactic WD can 
be smaller than the plotted one and, in any case, is direction dependent, we conclude that in 
the frequency interval ~ 10~^ — 10~^ Hz none of GW backgrounds of Galactic origin should be 
higher than the LISA sensitivity in the frequency bins A/ = 3 x 10^^ Hz. If LISA detects a 
GW background in this frequency interval, it is expected to be of a primordial origin. Although 
the low-frequency part of the binary confusion limit is somewhat more model-dependent (it is 
determined by the actual number of binary stars in the Galaxy, their space distribution and 
details of binary evolution), the calculated noise drops below the LISA sensitivity curve at 
frequencies below ~ 10~^ Hz. This leaves open for the search for cosmological backgrounds of 
primordial origin some low-frequency portion of the LISA sensitivity window, in addition to the 
already discussed interval ~ 10~^-10~^ Hz. 

5.2 Gravitational wave noise from extra-Galactic binaries 

Simple estimates show that the isotropic extra-Galactic background is expected to be one and 
a half order of magnitude smaller than the sky-averaged GW noise from Galactic binaries (see 



114, 33, 36, 0). Consider a volume of space with radius r = 300 Mpc. This is a large 
volume, but one can still neglect effects of curvature and cosmological time-dependence of star 
formation rate. According to Eq. (p!^, event rate in this volume TZy is related with the Galactic 



event rate TZq as TZy = 3 x 10 7?.g- Formula (|29D , written for extra-Galactic sources, should 



now contain TZy instead of the galactic rate TZ, and r = 300Mpc instead of the characteristic 
Galactic distance r = lOkpc. Combining the numbers, one finds the relationship between the 



36 



extra-Galactic amplitude h^{f) and the Galactic amplitude h]y{f): 

hfFif) « 5 X 10-2/i^(/). (44) 

Thus, the noise amplitude from extra-Galactic binaries is expected to be a factor 20 smaller 
than the noise amplitude from binaries in our Galaxy. 

More sophisticated calculations take into account the somewhat larger star formation rate 
at larger red shifts []122| , 115|. According to these studies, unresolved extra-Galactic binaries can 



contribute up to 10% of the mean Galactic noise. This is still smaller than the projected LISA 
sensitivity and thus presents no danger of contamination. The contribution of unresolved extra- 



Galactic binaries is shown in Fig. 13 to the right of the frequency ~ 10~ Hz, where Galactic 
binaries become resolvable and their noise contribution sharply drops down. An additional 
distinctive feature of GW backgrounds from sources in distant galaxies is a certain angular 
anisotropy of the background caused by inhomogeneities in the distribution of galaxies over the 
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6 Relic Gravitational Waves and Their Detection 
6.1 Introduction 

The existence of relic gravitational waves is a consequence of quite general assumptions. Es- 
sentially, we rely only on the validity of general relativity and basic principles of quantum field 
theory. The strong variable gravitational field of the early Universe amplifies the inevitable zero- 
point quantum oscillations of the gravitational waves and produces a stochastic background of 
relic gravitational waves measurable today [124, |125| , The detection of relic gravitational 



waves is the only way to learn about the evolution of the very early Universe, up to the limits of 
the Planck era and the big bang. It is important to appreciate the fundamental and unavoidable 
nature of this mechanism. Other physical processes can also generate stochastic backgrounds 
of gravitational waves. But those processes either involve many additional hypotheses, which 
may turn out to be not true, or produce a gravitational wave background (like the one from 
binary stars in the Galaxy) which should be treated as an unwanted noise rather than a useful 
and interesting signal. The scientific importance of detecting relic gravitational waves has been 
stressed on several occasions (see, for example, [Q, ||, |6|). 

The central notion in the theory of relic gravitons is the phenomenon of super-adiabatic 
(parametric) amplification. The roots of this phenomenon are known in classical physics, and 
we will remind its basic features. As every wave-like process, gravitational waves are amenable 
to the concept of a harmonic oscillator. The fundamental equation for a free harmonic oscillator 
is 

q + uj'^q = 0, (45) 

where q can be a displacement of a mechanical pendulum or a time-dependent amplitude of a 
mode of the physical field. The energy of the oscillator can be changed by an external force or, 
alternatively, by a parametric influence, that is, when a parameter of the oscillator, for instance 
the length of a pendulum, is being changed. In the first case, the fundamental equation takes 
the form 

q + i^\ = f{t), (46) 



whereas in the second case Eq. (45) takes the form 

g + cu2(t)g = 0. (47) 
Equations (^) and (^^ are profoundly different, both, mathematically and physically. 
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Figure 6: Parametric amplification, a) variation of tlie lengtli of tlie pendulum, b) increased 
amplitude of oscillations. 

Let us concentrate on the parametric influence. We consider a pendulum of length L oscillat- 
ing in a constant gravitational field g. The unperturbed pendulum oscillates with the constant 
frequency oj = \fgjL. Fig. ^ illustrates the variation of the length of the pendulum L(f) by an 
external agent, shown by alternating arrows. Since L(f) varies, the frequency of the oscillator 
does also vary: = g/L[t). The variation in L{t) need not be periodic, but it cannot be too 
slow (i.e., adiabatic) for the result of the process to be significant. Otherwise, in the adiabatic 
regime of slow variations, the energy of the oscillator E and its frequency uj do change slowly, 
but E/ijo remains constant, so one can say that the "number of quanta" E/fiu in the oscillator 
remains fixed. In other words, for the creation of new "particles - excitations", the characteristic 
time of the variation should be comparable with the period of the oscillator and the adiabatic 
behaviour should be violated. After some duration of the appropriate parametric influence, the 
pendulum will oscillate at the original frequency, but will have a signiflcantly larger, than before, 
amplitude and energy. This is shown in Fig. |6|b. Obviously, the energy of the oscillator has been 
increased at the expense of the external agent (pump field). For simplicity, we have considered 
a familiar case, when the length of the pendulum varies, while the gravitational acceleration g 
remains constant. Variation in g represents a gravitational parametric influence that would be 
in an even closer analogy with what we study below. 

A classical oscillator must have a non-zero initial amplitude for the amplification mechanism 
to work. Otherwise, if the initial amplitude is zero, the final amplitude will also be zero. Indeed, 
imagine the pendulum to be strictly at rest, hanging straight down. No variation in its length will 
cause the pendulum to oscillate and gain energy. In contrast, a quantum oscillator does not need 
to be excited from the very beginning. The oscillator can be initially in its quantum-mechanical 
vacuum state. The inevitable zero-point quantum oscillations are associated with the vacuum 
state energy ^fiu). One can imagine a pendulum hanging straight down, but fluctuating with a 
tiny amplitude determined by the "half of the quantum in the mode" . In the classical picture, 
it is this tiny amplitude of quantum-mechanical origin that is being parametrically amplifled. 

The Schroedinger evolution of a quantum oscillator depends crucially on whether the oscil- 
lator is being excited parametrically or by a force. Consider the phase diagram {q,p), where q is 
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the displacement and p is the conjugate momentum. The vacuum state is described by a circle 
at the center of the phase-space (see Fig. The mean values of q and p are each zero, but their 
variances (i.e. the zero-point quantum fluctuations) are not zeroes and are equal to each other. 
The magnitudes of the variances are represented by the radius of the circle at the center. Under 
the action of a force, the vacuum state evolves into a coherent state. The mean values of p and q 
have increased, but the variances are still equal and are described by the circle of the same size 
as for the vacuum state. On the other hand, under a parametric influence, the vacuum state 
evolves into a squeezed vacuum state. (For a recent review of squeezed states see, for example, 
|126| ] and references there.) Its variances for the conjugate variables q and p are significantly 
unequal and are described by an ellipse. As a function of time, the ellipse rotates with respect 
to the origin of the (g, p) diagram, and the numerical values of the variances oscillate too. The 
mean numbers of quanta in the two states, one of which is coherent and the other a squeezed 
vacuum, can be equal (similar to the coherent and squeezed states shown in Fig. ^) but the sta- 
tistical properties of these states are significantly different. Among other things, the variance of 
the phase of the oscillator in a squeezed vacuum state is very small (hence the name, squeezed). 
Graphically, this is reflected in the fact that the ellipse is very thin, so that that the uncertainty 
in the angle between the horizontal axis and the orientation of the ellipse is very small. This 
highly elongated ellipse can be regarded as a portrait of the gravitational wave quantum state 
that is being inevitably generated by parametric amplification, and which we will be dealing 
with below. 

A wave-field is not a single oscillator, it depends on spatial coordinates and time, and may 
have several independent components (polarization states). However, the field can be decom- 
posed into a set of spatial Fourier harmonics. In this way we represent the gravitational wave 
field as a collection of many modes, many oscillators. Because of the nonlinear character of 
the Einstein equations, each of these oscillators is coupled to the variable gravitational field of 
the surrounding Universe. For sufficiently short gravitational waves of experimental interest. 
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this coupling was especially effective in the early Universe, when the condition of the adiabatic 
behaviour of the oscillator was violated. It is this homogeneous and isotropic gravitational field 
of all the matter in the early Universe that played the role of an external agent - the pump field. 
The variable pump field acts parametrically on the gravity-wave oscillators and drives them into 
multi-particle states. Concretely, the initial vacuum state of each pair of waves with oppositely 
directed momenta evolves into a highly correlated state known as the two-mode squeezed vac- 



uum state [127, 128, |129| , |130f| . The strength and duration of the effective coupling depends on 
the oscillator's frequency. They all start in the vacuum state but get excited by various amounts. 
As a result, a broad spectrum of relic gravitational waves is being formed. This spectrum is 
accessible to our observations today. 

6.2 Cosmological Gravitational Waves 

In the framework of general relativity, a homogeneous isotropic gravitational field is described 
by the line element 

ds^ = c^dt^ - a\t)6ijdx'dxK 
It is more convenient to introduce a new time coordinate t] and to write ds^ in the form 

ds2 = (r]) [d??^ - Sijdx'dx^] . (48) 

In cosmology, the function a{t) (or a(??)) is called scale factor. In our discussion, it will represent 
the gravitational pump field. 

Cosmological gravitational waves are small corrections hij to the metric tensor. They are 
defined by the expression 

ds^ = a2(r?)[dr/2 - {5ij + hij)dx'dx^]. (49) 

The functions hij{ri,x.) can be expanded over spatial Fourier harmonics e*"'' and e~*"^, where 
n is a constant wave vector. In this way, we reduce the dynamical problem to the evolution 
of time-dependent amplitudes for each mode n. Among six functions hij there are only two 
independent (polarization) components. This decomposition can be made, both, for classical 
and for quantized field hy. In the quantum version, the functions hij are treated as quantum- 
mechanical operators. We will use the Heisenberg picture, in which the time evolution is carried 
out by the operators while the quantum state is fixed. This picture is fully equivalent to the 
Schroedinger picture, discussed in the Introduction, in which the vacuum state evolves into a 
squeezed vacuum state while the operators are time-independent. 

The Heisenberg operator for the quantized real field hij can be written as 



(50) 



where C is a constant which will be discussed later. The creation and annihilation operators 
satisfy the conditions [Cn, c^] = Ss's^'^i 

n — m), CjijO) — 0, where |0) (for each n and 5) is the 
fixed initial vacuum state discussed below. The wave number n is related with the wave vector 
n by n = (dijU^n^)^^"^ . The two polarization tensors Pij{n) (s = 1,2) obey the conditions 

PijU^ = 0, Pij6''^ = 0, PijP = 26ss', Pij{-n) = Pij{n). 
The time evolution, one and the same for all n belonging to a given n, is represented by the 

s 

complex time-dependent function hnif])- This evolution is dictated by the Einstein equations. 
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The nonlinear nature of the Einstein equations leads to the coupling of hniv) with the pump 

s 

field a(r/). For every wave number n and each polarization component s, the functions hniv) 
have the form 

hniv) = -T^iuniv) + Kiv)], (51) 

where Un{ri) and Vn{r]) can be expressed in terms of three real functions (the polarization index 
s is omitted): r„ - the squeeze parameter, (f>n - the squeeze angle, On ~ the rotation angle, 

Un = e*^" cosh r„, Vn = e"*^^""^'^") sinh r„. (52) 

The dynamical equations for Unii]) and Vn{rj) 

■ dun , .a' ^ .dvn , .a' ^ 

I— — =nun+i—v„, i——=nvn + i—u„ (53) 

arj a dry a 

lead to the dynamical equations governing the functions rn{rj), 4>n{il) and 0n(??) il29| , |130| ]: 

r'^ = — cos 2(/)„, cj)'j^ = —n sin 20„ coth 2r„, 9'^ = —n sin 2(j)n tanh r„, (54) 

a a a 

where ' = d/dr/, and the evolution begins from r„ = 0. This value of r„ characterizes the 
initial vacuum state |0) which is defined long before the interaction with the pump field became 
effective, that is, long before the coupling term a' /a became comparable with n. The constant C 
should be taken as C = \/167r Ipi where Ipi = (Gh/c^)^^'^ is the Planck length. This particular 
value of the constant C guarantees the correct quantum normalization of the field: energy ^huj 
per each mode in the initial vacuum state. The dynamical equations and their solutions are 
identical for both polarization components s. 

Equations (|53| ) can be translated into the more familiar form of the second-order differential 

equation for the function fJ-n{v) = '^niv) + ^n('?) = (^{v)hn{v) |124| , |125| , |^: 



0. (55) 



Clearly, this is the equation for a parametrically disturbed oscillator (compare with Eq. 
^)). In absence of the gravitational parametric influence represented by the term a"/a, the 
frequency of the oscillator defined in terms of 77-time would be a constant: n. Whenever the 
term a" /a can be neglected, the general solution to Eq. ( ^5|) has the usual oscillatory form 



finiv) = ^ne-™" + i?„e^"^ (56) 

where the constants An, Bn are determined by the initial conditions. On the other hand, 
whenever the term a" /a is dominant, the general solution to Eq. (p5|) has the form 

/ d?] 

/"n(r/) = Cna + Dna / (57) 

In fact, this approximate solution is valid as long as n is small in comparison with |a'/o|. This 
is more clearly seen from the equivalent form of Eq. ( p5| ) written in terms of the function hn{'r]) 

|]: 

< + + n^hn = 0. (58) 

a 
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Figure 8: Effective potential U{r]). 

For growing functions a(r/), tliat is, in expanding universes, the second term in Eq.(|57|) is usually 
smaller than the first one (see below), so that, as long as n <C a' /a, the dominant solution is the 
growing function finiv) = Cna{r]), and 

hn = const. (59) 

Equation (|55| ) can also be treated as the Schroedinger equation for a particle moving in the 
presence of an effective potential U{r]) = a" /a. In those situations that are normally considered, 
the potential U{ri) has a bell-like shape and forms a barrier (see Fig. |8|). When a given mode n is 
outside the barrier, its amplitude /i„ is adiabatically decreasing with time: oc e^^"'^ / a{ri) . This 
is shown in Fig. |8| by oscillating lines with decreasing amplitudes of oscillations. The modes with 
sufficiently high frequencies do not interact with, and stay above, the barrier. Their amplitudes 
hn behave adiabatically all the time. For these high-frequency modes, the initial vacuum state 
(in the Schroedinger picture) remains the vacuum state forever. On the other hand, the modes 
that interact with the barrier are subject to the super-adiabatic amplification. Under the barrier 
and as long as n < a' /a, the function /i„ stays constant instead of the decreasing adiabatically. 
For these modes, the initial vacuum state evolves into a squeezed vacuum state. 
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After having formulated the initial conditions, the present day behaviour of r„, On (or, 
equivalently, the present day behaviour of hn) is essentially all we need to find. The mean number 
of particles in a two- mode squeezed state is 2sinh^r„ for each s. This number determines the 
mean square amplitude of the gravitational wave field. The time behaviour of the squeeze angle 
(j)n determines the time dependence of the correlation functions of the field. The amplification 
(that is, the growth of r„) governed by Eq. (^) is different for different wave numbers n. 
Therefore, the present day results depend on the present day frequency v {u = cn/27ra) measured 
in Hz. 

In cosmology, the function H = a/ a = ca' /a? is the time-dependent Hubble parameter. The 
function I = c/H is the time-dependent Hubble radius. The time-dependent wavelength of the 
mode n is A = 27ra/n. The wavelength A has this universal definition in all regimes. In contrast, 
the V defined asv = cn/ 2'Ka has the usual meaning of a frequency of an oscillating process only in 
the short- wavelength (high-frequency) regime of the mode n, that is, in the regime where A <C /. 
As we have seen above, the qualitative behaviour of the solutions to Eqs. (|55|), (^), depends 
crucially on the comparative values of n and a'/a, or, in other words, on the comparative values 
of A(r/) and This relationship is also crucial for the solutions to Eq. ([5^), as we shall see 
now. 

In the short-wavelength regime, that is, during intervals of time when the wavelength A(??) 
is shorter than the Hubble radius l{ifj = a? /a' , the term containing n in (^^ is dominant. The 
functions 0n(??) and Onijj) are, (pn — — nirj -\- Tjn)^ On — 0n 

where ijn is a constant. The factor 
cos 2(j)n is a rapidly oscillating function of time, so the squeeze parameter r„ stays practically a 
constant. This is the adiabatic regime for a given mode. 

In the opposite, long-wavelength regime, the term n can be neglected. The function (pn is 
tan (77) const/a^ (ry) , and the squeeze angle quickly approaches one of the two values: 0„ = 
or = vr (analog of "phase bifurcation" [ [L31| , |132[ |). When the long- wavelength regime, for a 
given n, begins the squeeze parameter r„(ry) grows with time according to 



r„(r/) PS In 



a{r]) 



where is the value of a(r/) at r/^,. The final value of r„ is 



In 



a* 



(60) 



(61) 



where a** is the value of a{rj) at rj^^,, when the long- wavelength regime and amplification come to 
the end. It is important to emphasize that it is not a "sudden transition" from one cosmological 
era to another that is responsible for amplification, but the entire interval of the long-wavelength 
(non-adiabatic) regime. 

After the end of amplification, the accumulated (and typically large) squeeze parameter r„ 
stays approximately constant. The mode is again in the adiabatic regime. In course of the 
evolution, the complex functions Un{r)) + v^{r]) become practically real, and one has hniv) ~ 
hniv) ~ a^^" cospniv)- Every amplified mode n of the field (50) takes the form of a product 



of a function of time and a (random, operator-valued) function of spatial coordinates; the mode 
acquires a standing- wave pattern. The periodic dependence cos </>„(??) will be further discussed 
below. 

It is clearly seen from the fundamental equations (54), (55), ( |5^ ) that the final results depend 



only on a(r/). Equations do not ask us the names of our favorite cosmological prejudices, they ask 
us about the pump field a(r/). Conversely, from the measured relic gravitational waves, we can 
deduce the behaviour of a{rj), which is essentially the purpose of detecting the relic gravitons. 
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6.3 Cosmological Pump Field 

With the chosen initial conditions, the final numerical results for relic gravitational waves depend 
on the concrete behaviour of the pump field represented by the cosmological scale factor a{rj). 
We know a great deal about a(r/). We know that a{rj) behaves as 0(77) oc r/^ at the present 
matter-dominated stage. We know that this stage was preceded by the radiation-dominated 
stage in which a{rj) oc rj. At these two stages of evolution the functions a{rj) are simple power- 
law functions of rj. What we do not know is the function a(rj) describing the initial stage of 
expansion of the very early Universe, that is, before the era of primordial nucleo-synthesis. It is 
convenient to parameterize 0(77) at this initial stage also by power-law functions of 77. First, this 
is a sufficiently broad class of functions, which, in addition, allows us to find exact solutions to 
our fundamental equations. Second, it is known |124, |125| , |9| that the pump fields o(ry) which 



have power-law dependence in terms of r], produce gravitational waves with simple power-law 
spectra in terms of v. These spectra are easy to analyze and discuss in the context of detection. 

We model cosmological expansion by several successive eras. Concretely, we take a{r]) at the 
initial stage of expansion (i-stage) as 

a{v) = lo\v\'''^ (62) 

where grows from —00, and 1 + (3 < 0. We will show later how the available observational 
data constrain the parameters lo and (3. The z-stage lasts up to a certain t] = 7]i, rji < 0. To 
make our analysis more general, we assume that the i-stage was followed by some interval of 
the z-stage {z from Zeldovich). It is known that an interval of evolution governed by the most 
"stiff' matter (effective equation of state p = e) advocated by Zeldovich, leads to a relative 
increase of gravitational wave amplitudes |124| , |125| , It is also known that the requirement 
of consistency of the graviton production with the observational restrictions does not allow the 



'stiff' matter interval to be too long |124, 125, ^, 133 |. However, we want to investigate any 



interval of cosmological evolution that can be consistently included. In fact, the z-stage of 
expansion that we include is quite general. It can be governed by a "stiffer than radiation" 
1 134] matter, as well as by a "softer than radiation" matter. It can also be simply a part of the 



radiation-dominated era. Concretely, we take a(r/) at the interval of time from rji to some rjg 
(z-stage) in the form 

a{ri) = loa,{ri - rip)^+'^% (63) 

where l+/3s > 0. For the particular choice (3s = 0, the z-stage reduces to an interval of expansion 
governed by the radiation-dominated matter. Starting from rjg and up to 7]2 the Universe was 
governed by the radiation-dominated matter (e-stage). So, in this interval of evolution, we take 
the scale factor in the form 

a{r]) = loae{r] - 7]e). (64) 
And, finally, from f] = i]2 the expansion switched to the matter-dominated era (m-stage): 

a(r?) = loamiv - rjmf. (65) 



A link between the arbitrary constants participating in Eqs. ( |6^ ) - ( pq ) is provided by the 
conditions of continuous joining of the functions a{rj) and a'{i]) at points of transitions 771, rjs, 

m- 

We denote the present time by rjn (R from reception). This time is defined by the observation- 
ally known value of the present-day Hubble parameter H{r]R) and Hubble radius Ih = c/Hirjn). 
For numerical estimates we will be using lu k, 2x 10^^ cm. It is convenient to choose rjfi—rjm = 1, 
so that a(r]R) = 21h- The ratio 

a{riR)/a{m) = C2 
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is believed to be around C2 = 10'*. We also denote 

a{r]2)/a{r]s) = Cs , a{r]s) / a{r]i) = Ci 



With these definitions, all the constants participating in Eqs. ( |6^ ) - ( |65D (except parameters (3 
and fig which should be chosen from other considerations) are being expressed in terms of Ih^ 
C2, Cs) S'lid Ci- For example, 

I I + 
The important constant Iq is expressed as 

lo = blHC2' CfCi^, (66) 

where b = 2^^^ /\1 + Note that 6=1 for fj = —2. (This expression for may help 

to relate formulas written here with the equivalent treatment ]135| ] which was given in slightly 
different notations.) The sketch of the entire evolution a{rj) is given in Fig. |9|. 

We work with the spatially-flat models (|4^). At every instant of time, the energy density 
e(r/) of matter driving the evolution is related with the Hubble radius l{r]) by 

3 

'^^(^) = ]2(^' (67) 

where k = SvrG/c^. For the case of power-law scale factors 0(7/) oc rj^^^, the effective matter 
pressure p{r]) is related to the energy density e{ri) by the effective equation of state 

P = ] ~ ^n. (68) 
^3(1+/?) ^ ^ 

For instance, p = for /3 = 1, p = |e for /? = 0, p = — e for /3 = —2, and so on. Each interval of 
the evolution (|6^)-(|65|) is governed by one of these equations of state. 

In principle, the function a(r/) could be even more complicated than the one that we consider. 
It could even include an interval of early contraction, instead of expansion, leading to a "bounce" 
of the scale factor. In the case of a decreasing a{r\) the gravitational-wave equation can still be 
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Figure 10: Function a' /a for the scale factor from Fig. ^ 



analyzed and the amplification is still effective |124, 125, ^. However, the Einstein equations 
for spatially-flat models do not permit a regular "bounce" of o(ry) (unless e vanishes at the 
moment of "bounce"). Possibly, a "bounce" solution can be realized in alternative theories. 



such as string-motivated cosmologies |13(:, 137, 138]. For a recent discussion of spectral slopes 
of gravitational waves produced in "bounce" cosmologies, see [|139||. 



6.4 Solving Gravitational Wave Equations 

The evolution of the scale factor 0(77) given by Eqs. ( |6^ ) - ( |65[ ) and sketched in Fig. ^ allows us to 



calculate the function a' /a. This function is sketched in Fig. 10. In all theoretical generality, the 
left-hand-side of the barrier in Fig. |^ could also consist of several pieces, but we do not consider 
this possibility here. The graph also shows the important wave numbers nn, n2, ng, ni: nn 
marks the wave whose today's wavelength \{r]R) = 2Tia{r]ii) /nu is equal to the Hubble radius 
Ih today. With our parameterization a{i]ji) = 21h, this wave-number is tlh = 4tt. n2 marks the 
wave whose wavelength A(r/2) = 2Tra{r]2)/n2 at ry = 772 is equal to the Hubble radius 1{t]2) at t] = 
r]2. Since X{vR)/Km) = ("-2/?^^^)[a(^?R)/a(r72)] and l{r]R)/l{r]2) = [a{r]R)/a{r]2)][a{r]R)/a{r]2)]^/'^, 



this gives us n2/nH 
find 



1/2 



^2 _A 
— '32 ' 



Working out in a similar fashion other ratios, we 

(69) 



n2 Us 



Solutions to the gravitational wave equations exist for any a(r/). At intervals of power-law 



dependence a(r/), solutions to Eq. (55) have simple form of the Bessel functions. We could have 
found piece-wise exact solutions to Eq. ( |55[ ) and join them in the transition points. However, 
we will use a much simpler treatment, which is sufficient for our purposes. We know that the 
squeeze parameter r„ stays constant in the short-wavelength regimes and grows according to 
Eq. ( |60|) in the long-wavelength regime. All modes start in the vacuum state, that is, r„ = 
initially. After the end of amplification, the accumulated value ( |6l|) stays constant up to today. 
To find today's value of e"^" we need to calculate the ratio a^^{n)/a^{n). For every given n, 
the quantity a,, is determined by the condition A(ry*) = /(r/*), whereas a** is determined by the 
condition A(r/*^,) = l{r]^,^). 

Let us start from the mode n = ni. For this wave number we have a^, = a*^, = a(f/i), and 



therefore, 



0. The higher frequency modes, i.e. n > ni (above the barrier in Fig. QG), have 
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never been in the amphfying regime, so we can write 

e*"" = 1, n>ni. (70) 

Let us now consider the modes n in the interval rii > n > Us- For a given n we need to know 
a*(n) and a^,^,{n). Using Eq. (|6^ one has a^,{n)/a^,{nl) = (ni/n)^^^ , and using Eq. ( |63|) one 
finds a^^:{n)/a^:^:{ns) = {ng/nY^I^" . Therefore, one finds 



n 



a*(n) a*(ni) \ n J \ni , 
Since a**(ns) = a{r]s), a*(ni) = a(r/i), and a{r]s) / a{rii) = Ci = (f^i/iT's)^^^'' , we arrive at 

a*{n) \ni) 

Repeating this analysis for other intervals of the decreasing n, we come to the conclusion that 



V 



n 

niJ 



ni > n > Us, 



IT'S > n > n2, 



e- = I^Y-'(!^X(!^) ,n2>n>n,. (71) 



UsJ \niJ 
n2j \niJ \ni, 

The mnemonic rule of constructing e*"" at successive intervals of decreasing n is simple. If the 
interval begins at Ux, one takes {n/nx)^*~^** and multiples with e*""^, that is, with the previous 
interval's value of e*"" calculated at the end of that interval n^- For the function a' /a that we 
are working with, the is always /?, whereas the takes the values j3s, 0, 1, at the successive 
intervals. 

The modes with n < nn are still in the long-wavelength regime. For these modes, we should 
take a{rip>) instead of a**(n). Combining with a*(n), we find 

M^"V!l£)'-VM''fM -'■,„<„„, (72) 

nnJ \n2j \niJ \niJ 

Formulas (|7^) - (72) give approximate values of r„ for all n. The factor e''" obeys the following 
inequalities is e''" > 1 for n < ni, and e^" S> 1 for n <^ rii, and determines the mean square 
amplitude of the gravitational waves. 

The mean value of the field hij is zero at every moment of time rj and in every spatial point 
x: (0|/iij(?7, x)|0) =0. The variance 

(0|/ii,(r/,x)/i^^(r?,x)|0) ^ {h^) 

is not zero, and it determines the mean square amplitude of the generated field - the quantity 
of interest for the experiment. Taking the product of two expressions (5C) one can show that 

oo 9 OO 

{h') = ^ J nY.\hn{v)\ dn^Jh\n,v)—. (73) 

n s=l n 



Using the representation (|5l|), (52) in Eq. ( [73|) one can also write 



C f 

Ut^) = „ „ / ndn(cosh2r„ + cos2(/)„sinh2r„). (74) 
TT'^a^ J 
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We can now consider the present era and use the fact that e"^" are large numbers for ah n in the 
interval of our interest ni > n > uh- Then, we can derive 

h{n,r]) — / . rae''" cos(/>rt(f?) = 8^ fr^) ( —] cos . (75) 

TT a{r]ii) \Ih J KuhJ 

The quantity h{n, rf) is the dimensionless spectral amplitude of the field whose numerical value 
is determined by the calculated squeeze parameter r„. The oscillatory factor cos(/>„(r/) reflects 
the squeezing (standing wave pattern) acquired by modes with ni > n > uh. For modes 
with n < UH this factor is approximately 1. For high-frequency modes n ^ nn one has 
(f>n{'n) ~ ^(^ ~ Vn) 2> 1) so that h{n,rj) makes many oscillations while the scale factor a{rj) is 
practically fixed at a{rjii). 

The integral (|7j) extends formally from to oo. Since r„, ~ for n > ni, the integral diverges 
at the upper limit. This is a typical ultra-violet divergence. It should be discarded (renormalized 
to zero) because it comes from the modes which have always been in their vacuum state. At the 
lower limit, the integral diverges, if /? < —2. This is an infra-red divergence which comes from 
the assumption that the amplification process has started from infinitely remote time in the past. 
One can deal with this divergence either by introducing a lower frequency cut-off (equivalent to 
the finite duration of the amplification) or by considering only the parameters f3 > —2, in which 
case the integral is convergent at the lower limit. It appears that the available observational 
data (see below) favour this second option. The particular case /? = —2 corresponds to the de 
Sitter evolution 0(77) oc |?7|~^. In this case, the h{n) found in Eqs. (^), (fz^) does not depend 
on n. This is known as the Harrison-Zeldvich, or scale-invariant, spectrum. 

The spectral amplitudes h{n) can also be derived using the approximate solutions (|^), (|57|) 
to the wave equation (^5|). This method gives exactly the same, as in Eqs. ([75|), (|70|) - ([7^) 
numerical values of h{n), but does not reproduce the oscillatory factor cos(j)niv)- 

One begins with the initial spectral amplitude hi{n) defined by quantum normalization: 
hi{n) = 8^/^T{lpl/Xi). This is the amplitude of the mode n at the moment r]^ of entering the 
long wavelength regime, i.e. when the mode's wavelength Aj is equal to the Hubble radius /(r/*). 
For Aj one derives 

Thus, we have 

y ^2+|3 

h,{n) = A( — ) , (77) 

where A denotes the constant 

A = 68V^^ . (78) 
The numbers hi{n) are defined at the beginning of the long- wavelength regime. In other words. 



they are given along the left-hand-side slope of the barrier in Fig. 10. We want to know the final 
numbers (spectral amplitudes) h{n) which describe the field today, at r]ji. 

According to the dominant solution hn{i]) = const, of the long- wavelength regime (see Eq. 
(|5^), the initial amplitude hi{n) stays practically constant up to the end of the long- wavelength 
regime at Tj^^^ that is, up to the right-hand-side slope of the barrier. (The second term in 
Eq. ( ^7|) could be important only at the z-stage and only for parameters [5s < —(1/2), which 
correspond to the effective equations of state p > e. In order to keep the analysis simple, we do 
not consider those cases.) After the completion of the long- wavelength regime, the amplitudes 
decrease adiabatically in proportion to l/a{rj), up to the present time. Thus, we have 

,^„>^A(^r''-r^. ,79) 
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Let us start from the lower end of the spectrum, n < nn, and go upward in n. The modes 
n < uh have not started yet the adiabatic decrease of the amphtude, so we have 



h(n) = A{ — , n<nH- (80) 

Now consider the interval n2 > n > nn- At this interval, the a^*{n) / airm) scales as [nulnf'-, 
so we have ^ 

h(n) = A\—\ , n2>n>nH. (81) 
\nHj 

At the interval Ug > n > n2 the ratio a<t:^{n) / a{r]R) = [a**(n)/a(r/2)][a(r/2)/a(?7_R)] scales as 
{n2/n){nH /n2)'^ , so we have 

h{n)=A[ — — , ns>n>n2. (82) 
KuhJ n2 

Repeating the same analysis for the interval ni > n > we find 



h{n)=A\ — — ni>n>ns. (83) 

ynuJ KuhJ n2 

It is seen from Eq. (^) that an interval of the z-stage with /3s < (the already imposed 
restrictions require also (—1/2) < (3s) bends the spectrum h{n) upwards, as compared with Eq. 
(|8^), for larger n. If one recalls the relationship ( |66|) between and Ijj and uses ([n]), ([7^ ) in 
Eq. ( [75| ) one arrives exactly at Eqs. (|80D -(p^ ) up to the oscillating factor cos(/>„(r/). 

Different parts of the barrier in Fig. |lO| are responsible for amplitudes and spectral slopes 
at different intervals of n. The sketch of the generated spectrum h{n) in conjunction with the 



form of the barrier is shown in Fig. |ll 



The present day frequency of the oscillating modes, measured in Hz, is defined as = 
cn/2TTa{r]ji). The lowest frequency (Hubble frequency) is i^h = c/Ih- For numerical estimates 
we will be using Vfj ~ lO"^*^ Hz. The ratios of n are equal to the ratios of u, so that, for example, 
n/nu = V jvu. For high-frequency modes we will now often use the ratios of v instead of ratios 
of n. 

In addition to the spectral amplitudes h(n) the generated field can be also characterized by 
the spectral energy density parameter Q.g(n). The energy density of the gravitational wave 
field is ^ ^ 

h^^nhan = —rh^^ ha'. 



^^9 ~ /I fl'^ijfi 



4a2 



The mean value (0|e„(r/, x)|0) is given by 



2 



1 C [ \s' 







For high-frequency modes, it is only the factor e^*"'' that needs to be differentiated by r]. After 

2 



averaging out the oscillating factors, one gets 



s' 
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K 


= n 


hn 



so that 



oo 

K{e,) = ^ / n^h\n)^. (85) 


In fact, the high-frequency approximation, that has been used, permits integration over lower n 
only up to uh- And the upper limit, as was discussed above, is in practice ni, not infinity. The 
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Figure 11: Amplitudes and spectral slopes of h{n) are determined by different parts of the 
barrier a' /a. 
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parameter Qg is defined as Qg = {eg)/e, where e is given by Eq. ( p7\j (critical density). So, we 
derive 



and 

^,{y) = '^h\v){^^y . (86) 

The dimensionless quantity ^gif) is useful because it allows us to quickly evaluate the 
cosmological importance of the generated field in a given frequency interval. However, the 
primary and more universal concept is h{v), not ^g{v). It is the field, not its energy density, 
that is directly measured by the gravity-wave detector. One should also note that some authors 
use quite a misleading definition 0,g{f) = {1/ pc){dpgw/dln f) which suggests differentiation 
of the gravity-wave energy density by frequency. This would be incorrect and could cause 



disagreements in numerical values of Jig. Whenever we use Qg(i'), we mean relationship (8£); 



and for order of magnitude estimates one can use [124, 125, 



f),(^)«/.2(^)(^-^)'. (87) 
6.5 Theoretical and Observational Constraints 

The entire theoretical approach is based on the assumption that a weak quantized gravity-wave 
field interacts with a classical pump field. We should follow the validity of this approximation 
throughout the analysis. The pump field can be treated as a classical gravitational field as long 
as the driving energy density e is smaller than the Planck energy density, or, in other words, as 
long as the Hubble radius l{rj) is greater than the Planck length Ipi. This is a restriction on the 
pump field, but it can be used as a restriction on the wavelength Aj of the gravity-wave mode n 
at the time when it enters the long- wavelength regime. If /(?]*) > Ipi, then Aj > Ipi. The Aj is 
given by Eq. ([76|). So, we need to ensure that 



At the lowest-frequency end u = vu this inequality gives b{lpi/lo) < 1. In fact, the observational 
constraints (see below) give a stronger restriction: 

6^ « lo-^ 

which we accept. Then, at the highest-frequency end u = ui we need to satisfy 

Let us now turn to the generated spectral amplitudes h^v). According to Eq. ( |80[ ) we have 
h{vH) ~ b8y/TT{lpi/lo). The measured microwave background anisotropics, which we discuss 
below, require this number to be at the level of 10~^, which gives the already mentioned Eq. 
(|88|). The quantity /i(i^i) at the highest frequency ui is given by Eq. (^): 



< 10^ (89) 



Vh) V2 ' 
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Using Eq. ( |66D this expression for h{ui) can be rewritten as 

%i) = 8V^^^ = 8V^^, (90) 

where Ai = c/ui. This last expression for /i(z^i) is not surprising: the modes with v '>vi are stiU 
in the vacuum state, so the numerical value of /i(i^i) is determined by quantum normalization. 

All the amplified modes have started with small initial amplitudes hi, at the level of zero- 
point quantum fluctuations. These amplitudes are also small today, since the hi could only stay 
constant or decrease. However, even these relatively small amplitudes should obey observational 
constraints. We do not want the in the high-frequency modes, which might affect the rate 
of the primordial nucleosynthesis, to exceed the level of 10^^. This means that Q.g{vi) cannot 
exceed the level of 10~^ or so. The use of Eq. ( |86| ) in combination with Qg{vi) ~ 10~^ and 
h{vi) from Eq. (|90|), gives us the highest allowed frequency i/i ~ 3 x 10^'^ Hz. We will use this 
value of vi in our numerical estimates. Returning with this value of ui to Eq. ( |89| ) we find that 
the parameter (3 can only be (3 < —1.8. We will be treating (3 = —1.8 as the upper limit for the 
allowed values of [3. 

We can now check whether the accepted parameters leave room for the postulated z-stage 
with 13s < 0. Using Eq. ( |66| ) we can rewrite Eq. ( [8^ ) in the form 

fil)-7fl)^-il. (91) 

Ipi \T^H/ \^s/ 

We know that U2/i^h = 10^ and vxjvs is not smaller than 1. Substituting all the numbers in Eq. 



(91) one can find that this equation cannot be satisfied for the largest possible /3 = —1.8. In 
the case /? = —1.9, Eq. ( pi] ) is only marginally satisfied, in the sense that a significant deviation 
from /3s = toward negative /3s can only last for a relatively short time. For instance, one can 
accommodate /J^ = —0.4 and = 10^ Hz. On the other hand, if one takes /3 = —2, a somewhat 
longer interval of the z-stage with /3s < can be included. For instance, Eq. (|9l|) is satisfied if 
one accepts Ug = 10^^ Hz and f3s = —0.3. This allows us to slightly increase h{iy) in the interval 
Us < 1^ < 1^1, SiS compared with the values of /i(z^) reached in the more traditional case /? = —2, 
/3s = 0. In what follows, we will consider consequences of this assumption for the prospects 
of detection of the produced gravitational wave signal. Finally, let us see what the available 
information on the microwave background anisotropies |140| , 141 1 allows us to conclude about 
the parameters /? and lo- 

Usually, cosmologists operate with the spectral index n (not to be confused with the wave 
number n) of primordial cosmological perturbations. Taking into account the way in which 
the spectral index n is defined, one can relate n with the spectral index /3 -|- 2 that shows 
up in Eq. (|80|). The relationship between them is n = 2(3 + 5. This relationship is valid 
independently of the nature of cosmological perturbations. In particular, it is valid for density 
perturbations, in which case the h{n) of Eq. ( |80|) is the dimensionless spectral amplitude of 
metric perturbations associated with density perturbations. If primordial gravitational waves 
and density perturbations were generated by the mechanism that we discuss here (an assumption 
that is likely to be true) then the parameter /? that participates in the spectral index is the same 
as the one that participates in the scale factor of Eq. (^) . Primordial gravitational waves and 
primordial density perturbations with the same spectral index P + 2 produce approximately the 
same lower-order multiple distributions of large-scale anisotropies. 

The evaluation of the spectral index n of primordial perturbations have resulted in n = 



1.2 lb 0.3 [141] or even in a somewhat higher value. A recent analysis [142| of all available data 
favors n = 1.2 and the quadrupole contribution of gravitational waves twice as large as that 
of density perturbations. One can interpret these evaluations as an indication that the true 
value of n lies somewhere near n = 1.2 (hopefully, the planned new observational missions will 
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determine this index more accurately). This gives us the parameter f3 somewhere near j3 = —1.9. 
We wiU be using f3 = —1.9 in our estimates below, as the observationally preferred value. The 
parameter (3 can be somewhat larger than /? = —1.9. However, as we already discussed, the 
value (3 = —1.8 (n = 1.4) is the largest one for which the entire approach is well posed. The 
Harrison-Zeldovich spectral index n = 1 corresponds to /3 = —2. 

The observed quadrupole anisotropy of the microwave background radiation is at the level 
6T/T « 10^^. The quadrupole anisotropy that would be produced by the spectrum ( |80| ) - ( |83| ) is 
mainly accounted for by the wave numbers near nn- Thus, the numerical value of the quadrupole 
anisotropy produced by relic gravitational waves is approximately equal to A. According to 
general physical considerations and detailed calculations [|143| , |144| , |145| ] , the metric amplitudes 
of long-wavelength gravitational waves and density perturbations generated by the discussed 
amplification mechanism are of the same order of magnitude. Therefore, they contribute roughly 
equally to the anisotropy at lower multipoles. This gives us the estimate A ~ 10~^, that we 



have already used in Eq. (p8D. It is not yet proven observationally that a significant part of 
the observed anisotropics at lower multipoles is indeed provided by relic gravitational waves, 
but we can at least assume this with some degree of confidence. It is likely that the future 
measurements of the microwave background radiation, and especially its polarization, will help 
us verify this theoretical conclusion. 

Combining all the evaluated parameters together, we show in Fig. 12 the expected spectrum 
of /i(z^) for the case /? = —1.9. A small allowed interval of the z-stage is also included. The 
intervals of the spectrum accessible to space- and ground-based interferometers are indicated by 
vertical lines. 



It is necessary to note |143, 144, 145, |146| | that a confirmation of any n > 1 (/? > —2) 
will mean that the very early Universe was not driven by a scalar field - the cornerstone of 
inflationary considerations. Indeed, the n > 1 (/3 > —2) requires the effective equation of state 
at the initial stage of expansion to be e + p < (see Eq. (|68|) ). But this requirement cannot be 
accommodated by any scalar fleld ip{t) whatever the scalar field potential V{(p) may be. The 
energy density of the scalar field is e = 0^/2 + V{(p), whereas its pressure is p = — V{ip), 
so that e + p cannot be negative. The available data do not prove yet that n > 1, but this 
possibility seems likely. 

It is also necessary to comment on a certain damage to gravitational wave research that 
was inflicted by the so called "standard inflationary result" . The "standard inflationary result" 
predicts infinitely large amplitudes of density perturbations 5p/p in the interval of spectrum 
with the Harrison-Zeldovich slope n = 1 (/3 = —2): 6p/p oc V^^'^/V' oc — n. The metric 

(curvature) perturbations hs accompanying 6p/p are also predicted to be infinitely large, in the 
same proportion. Infiationary literature conceals this predicted infinity of density perturbations 
by writing the ratio of the gravitational wave amplitude hx to the scalar metric amplitude hs- 
hx/hs ~ 7\/l — n, and declaring that the contribution of gravitational waves to the cosmic 
microwave background (CMB) anisotropics should be zero, or almost zero, for n ~ 1. Thus, the 
"standard" inflationary theory shifts the spectrum of relic gravitational waves, similar in shape 



to the one shown in Fig. |12|, down by many orders of magnitude. This claim has incorrectly 
led to discarding the gravitational wave contribution in the analysis of CMB data. Although 
one of the recent best fits to the available data indicates the presence of the gravitational wave 
contribution [147] this fit has been ignored. For many years, infiationary theorists claimed 
that their arbitrarily large density perturbations were caused by the "big amplification during 
reheating". It is now universally accepted that this explanation is false. The scalar metric 
perturbation, similarly to the gravitational wave perturbation (see Eq. (|59|)), remains constant 
during the long wavelength regime, that is, its numerical value does not change on the way 
from the first "Hubble-radius-crossing" to the second "Hubble-radius-crossing" . This fact is 
also refiected in the constancy of the so-called "conserved" gauge- invariant quantity (, which is 
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Figure 12: Expected spectrum /i(z^) for the case /? = —1.9. 
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in the center of inflationary analysis (see, for example, |148| ]). Since the "standard" inflationary 
theory predicts arbitrarily large numerical value of and this quantity is "conserved" during the 
evolution, this arbitrarily large number must have been postulated from the very beginning, as 
quantum normalization. This happens if one incorrectly assigns the quantum normalization to 
the scalar field perturbations alone, as if they were a free scalar test field in the De-Sitter space- 
time. Then, the quantity C, calculated from the perturbed Einstein equations, is arbitrarily 
large at the first "Hubble-radius-crossing" , and this number is being transmitted to the second 
"Hubble-radius-crossing" . In reality, however, scalar field perturbations are always coupled to 
metric (gravitational field) perturbations. The correct quantization of the combined degree of 
freedom renders the scalar metric (curvature) perturbation hs finite and small, and of the same 
order of magnitude as gravitational wave perturbation hx- The inflationary formula for 6p/p 



is incorrect in that it misses the dimensionless factor y—H/H'^ which cancels out the zero in 
the denominator and makes the generated density perturbations finite and small, even in the 
interval of spectrum with the Harrison-Zeldovich slope. The "standard inflationary result" is in 
a severe conflict not only with theory but with observations too: when the observers marginalize 
their data to n = 1 (enforce this value of n in data analysis) they find finite and small density 
perturbations, instead of infinitely large perturbations predicted by inflationary theorists. [For 
analytical expressions of the "standard inflationary result" see inflationary articles, including 
recent reviews. For graphical illustration of the divergent density perturbations and quadrupole 
anisotropics, predicted by inflationary theorists, see, for example, |149| ]. For critical analysis and 



disagreement with the "standard inflationary result" see [143, 144, 145|]. The most recent articles 



dealing with perturbations in quasi-deSitter models rightly emphasize the expected substantial 



contribution of gravitational waves to the large scale CMB anisotropics [150|. In short, general 



relativity and quantum field theory do not produce the "standard inflationary result" . 
6.6 Detectability of Relic Gravitational Waves 

We switch now from cosmology to prospects of detecting the predicted relic gravitational waves. 
The ground-based [|^, 14, 151 1 and space-based [16, |152| ] laser interferometers (see also [ 153 , 



154, |155|[ ) will be in the focus of our attention. We use laboratory frequencies u and intervals of 
laboratory time t {cdt = o(ry/j)dr/). Formulas (|82|) and (^3|), with A = 10~^, = 10^, and 

the oscillating factor restored, can be written as 



h{v, t) K, IQ-'^ cos[2T:u{t - t^)] ( — j , V2<u<Vs (92) 



and 

h{v, t) w 10"^ cos[2'nv{t - ty)] — — ■ Us < v < vi (93) 

where the deterministic (not random) constant ty does not vary significantly from one frequency 
to another at the intervals Az/ ~ u. The explicit time dependence of the spectral variance 
h'^{v,t) of the field, or, in other words, the explicit time dependence of the (zero-lag) temporal 
correlation function of the field at every given frequency, demonstrates that we are dealing 
with a non-stationary process (a consequence of squeezing and severe reduction of the phase 
uncertainty) . We will first ignore the oscillating factor and will compare the predicted amplitudes 
with the sensitivity curves of advanced detectors. The potential reserve of improving the SNR 
by exploiting squeezing will be discussed later. 



Let us start from the Laser Interferometer Space Antenna (LISA) [16[. The instrument will 
be most sensitive in the interval, roughly, from 10^'^ Hz to 10^^ Hz, and will be reasonably 
sensitive in a broader range, up to frequencies 10~^ Hz and 1 Hz. The sensitivity graph of LISA 
to a stochastic background is usually plotted under the assumption of a 1-year observation time. 
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Figure 13: Expected spectrum /? = —1.9 and other possible spectra in comparison with the 
LISA sensitivity. 



that is, the root-mean-square (r.m.s.) instrumental noise is being evaluated in frequency bins 
Ai/ = 3 X 10^*^ Hz around each frequency v. We need to rescale our predicted amplitude h{y) 
to these bins. 



The mean square amphtude of the gravitational wave field is given by the integral (73). 
Thus, the r.m.s. amplitude in the band Ai' centered at a given frequency v is given by the 
expression 

h{u,Ai^) = h{u)J—. (94) 




We use Eqs. (|92D, (pSj ) and calculate expression ( p4\ j assuming Au = 3 x 10~ Hz. The results 
are plotted in Fig. ^ Formula (^) has been used throughout the covered frequency interval 
for the realistic case /? = —1.9 and for the extreme case /? = —1.8. The line marked z-model 
describes the signal produced in the composite model with /? = — 2 up to i^s = 10"'' Hz (formula 
(p^) and then followed by formula (^) with Ps = —0.3. This model gives the signal a factor of 
3 higher at u = Hz, than the model /? = —2 extrapolated down to this frequency. 

There is no doubt that the signal /? = —1.8 would be easily detectable even with a single 
instrument. The signal /3 = —1.9 is marginally detectable, with the SNR around 3 or so, in a 
quite narrow frequency interval near and above the frequency 3 x 10^^ Hz. However, at lower 
frequencies one would need to be concerned with the possible gravitational wave noise from 
unresolved binary stars in our Galaxy (see Section ^). The further improvement of the expected 
LISA sensitivity by a factor of 3 may prove to be crucial for a confident detection of the predicted 
signal with /3 = —1.9. 

Let us now turn to the ground-based interferometers operating in the interval from 10 Hz to 
10^ Hz. The best sensitivity is reached in the band around v = 10^ Hz. We take this frequency 
as the representative frequency for comparison with the predicted signal. We will work directly 
in terms of the dimensionless quantity h{i'). If necessary, the r.m.s. amplitude per Hzi/2 at a 
given 1/ can be found simply as h{y)l\fv. The instrumental noise will also be quoted in terms 
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Figure 14: Full spectrum /i(z^) accessible to laser interferometers. 



of the dimensionless quantity hex{^)- 

The expected sensitivity of the initial instruments at z/ = 10^ Hz is h^x = 10^^^ or better. The 
theoretical prediction at this frequency, following from (|9^), (|93|) with j3s = 0, is hth = 10^^^ for 
(5 = —1.8, and hth = 10"^^ for j3 = —1.9. Therefore, the gap between the signal and noise levels 
is from 2 to 4 orders of magnitude. The expected sensitivity of the advanced interferometers, 
such as LIGO-II |l^, can be as high as hex = 10~^^. In this case, the gap vanishes for the 



j3 = —1.8 signal and reduces to 2 orders of magnitude for the /? = —1.9 signal. Fig. 14 illustrates 
the expected signal in comparison with the LIGO-II sensitivity. Since the signal lines are plotted 
in terms of h{v)^ the LISA sensitivity curve (shown for periodic sources) should be raised and 
adjusted in accordance with Fig. 

A signal below noise can be detected if the outputs of two or more detectors can be cross 
correlated. (For the early estimates of detectability of relic gravitational waves see [156|.) The 
cross correlation will be possible for ground-based interferometers, several of which are currently 
under construction. The gap between the signal and the noise levels should be covered by a 
sufficiently long observation time r. The duration r depends on whether the signal has any 
temporal signature known in advance, or not. We start from the assumption that no temporal 
signatures are known in advance. In other words, we first ignore the squeezed nature of the 
relic background and work under the assumption that the squeezing cannot be exploited to our 
advantage. 

The response of an instrument to the incoming radiation is s{t) = Fijh^^ where Fij depends on 
the position and orientation of the instrument. Since the h^^ is a quantum- mechanical operator 
(see Eq. (|50|)) we need to calculate the mean value of a quadratic quantity. The mean value of 
the cross correlation of responses from two instruments (0|si(t)s2(i)|0) will involve the overlap 
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reduction function |157, 158, |159| , 160 1, which we assume to be not much smaller than 1 [159]. The 
SNR in the measurement of the amphtude of a signal with no specific known features increases 
as (ri/)^/^, where v is some characteristic central frequency (for more detail see Sec. ^). 

We apply the guaranteed law (jv)^!'^ to initial and advanced instruments at the represen- 
tative frequency v = 10^ Hz. This law requires a reasonably short time r = 10^ sec in order 
to improve the SjN in initial instruments by two orders of magnitude and to reach the level 
of the signal with extreme spectral index /? = —1.8. The longer integration time or a better 
sensitivity will make the SjN larger than 1. In the case of a realistic spectral index /3 = —1.9 
the remaining gap of 4 orders of magnitude can be covered by the combination of a significantly 
better sensitivity and a longer observation time (not necessarily in one non-interrupted run). 
The sensitivity of the advanced laser interferometers, such as LIGO 11, at the level h^.^ = 10~^^ 
and the same observation time r = 10^ sec would be sufficient for reaching the level of the 
predicted signal with /? = —1.9. 

An additional increase of S/N can be achieved if the statistical properties of the signal can 
be properly exploited. Squeezing is automatically present at all frequencies from uh to i^i. The 
squeeze parameter r is larger in gravitational waves of cosmological scales, and possibly the 



periodic structure in Eq. (75) can be better revealed at those scales. However, we are interested 
here in frequencies accessible to ground based interferometers, say, in the interval 30-100 Hz. If 
our intention were to monitor one given frequency z/ from the beginning of its oscillating regime 
and up till now, then, in order to avoid the destructive interference from neighbouring modes 
during all that time, the frequency resolution of the instrument should have been incredibly 
narrow, of the order of 10~^^ Hz. Certainly, this is not something what we can, or intend to 
do. Although the amplitudes of the waves have adiabatically decreased and their frequencies 
redshifted since the beginning of their oscillating regime, the general statistical properties of the 
discussed signal are essentially the same now as they were 10 years after the big bang or will be 
1 million years from now. 

The periodic structure (^) may survive at some level in the instrumental window of sensi- 
tivity from fmm (minimal frequency) to I'max (maximal frequency). The mean square value of 
the field in this window is 

I h^{v,t)— = lQ-^^—^ I cos^[2Tiu{t-t^)y'^+^dv . (95) 

J V Vli J 

Because of the strong dependence of the integrand on frequency, z^~^'^ or i/"^ *^, the value of the 
integral (|95| ) is determined by its lower limit. Apparently, the search through the data should 
be based on the periodic structure that may survive at v = Vmin- As an illustration, one can 
consider such a narrow interval /S.v = Vmax — '^min that the integral (|9^) can be approximated 
by the formula 



'-^max 



/ h {U,t) ^ 10 cos [2TTVmin[t - tmin)\ ■ 

J V \ J \VminJ 

Clearly, the correlation function is strictly periodic and its structure is known in advance, in 
contrast to other possible signals. This is a typical example of using the a priori information. 
Ideally, the gain in S/N can grow as (ri/mm)^^^- This would significantly reduce the required 
observation time r. For a larger Az^, even an intermediate gain between the guaranteed law 
(ri/)^/^ and the law (ri^)^/^, from the matched filtering technique, would help. This could 
potentially make the signal with /? = —1.9 measurable even by the initial laser interferometers. 
A straightforward application of (Psl) exploiting squeezing may not be possible, as argued in a 



recent study [|161| , but more sophisticated methods are not excluded. 
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For frequency intervals covered by bar detectors and electromagnetic detectors, the expected 
results follow from the same formulas (92), (93) and have been briefly discussed elsewhere 
1 15^ , [135| ]. 



6.7 Conclusion 

It would be strange, if the predicted signal at the level corresponding to [i = —1.9 were not seen 
by the instruments capable of its detection. There aren't so many cosmological assumptions 
involved in the derivation, that could prove wrong, thus invalidating our predictions. On the 
other hand, it would be even more strange (and even more interesting) if the relic gravitational 
waves were detected at the level above the (3 = —1.8 line. This would mean that there is 
something fundamentally wrong in our basic cosmological premises. In the most favorable case, 
the detection of relic gravitational waves can be achieved by the cross-correlation of outputs 
of the initial laser interferometers in LIGO, VIRGO, GEO600. In the more realistic case, the 
sensitivity of advanced ground-based and space-based laser interferometers will be needed. The 
specific statistical signature of relic gravitational waves, associated with the phenomenon of 
squeezing, is a potential reserve for further improvement of the SNR. 



7 Gravitational wave detectors and their sensitivity 
7.1 Current status of GW antennas 

Currently, there are a number of bar detectors in operation: some of these operate at room- 
temperature and some others at cryogenic temperatures. Bar detectors are resonant, narrow- 
band detectors. They can detect signal amplitudes h ~ lO"^*' in a band width of 10-20 Hz around 
a central frequency of 1 kHz. Asymmetric supernovae in our Galaxy are the best candidate 
sources for these detectors. They may also see continuous radiation emitted by a neutron star 
if the frequency happens to lie in their sensitivity band. 

Interferometric detectors currently under construction will increase our ability to directly 
observe gravitational waves. The Japanese have already built a 300 m detector in Tokyo, 



Japan |162]. Several other projects are now nearing completion: The British-German col- 



laboration is constructing a 600 m interferometer (GEO) in Hannover, Germany jl63|] , the 



French-Italian collaboration is building a 3 km detector (VIRGO) near Pisa, Italy [ 164j| and the 
Americans are building two 4 km antennas (LIGO), one in Livingston and the other in Hanford 
|l5|| in the U.S.A. These detectors will start taking data between 2001 and 2003. The larger 
of these detectors, LIGO and VIRGO, are likely to be upgraded in sensitivity by an order of 
magnitude with a better low- frequency performance in 2005. These ground based interferom- 
eters will eventually be sensitive to sources in the frequency range from 10 Hz to several kHz. 
In addition to these ground based antennas, there is a plan to place an interferometer in space 
by the end of this decade. The Laser Interferometer Space Antenna (LISA) consists of three 
drag- free satellites, forming an equilateral triangle of size 5 million km, in a heliocentric orbit, 
lagging behind the Earth's orbit by 20°. LISA will be sensitive to waves in the low- frequency 
band of lO'^^-lO'^ Hz. 



7.2 Sensitivity of a GW antenna 

The performance of a GW detector is characterized by the power spectral density (henceforth 
denoted PSD) of its noise background. One can construct the noise PSD as follows: A GW 
detector outputs a dimensionless data train, say x{t), which in the case of an interferometer is 
the relative strain in the two arms. In the absence of any GW signal the detector output is just 
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an instance of noise n(t), that is, x{t) = n{t). The noise auto-correlation function n is defined 

asB 

K{ti,t2)=n{ti)n{t2), (96) 

where an overhne indicates the average over an ensemble of noise realizations. In general, k 
depends both on ti and t2- However, if the detector output is a stationary noise process, i.e. its 
performance is, statistically speaking, independent of time, then k depends only on r = t2 — ^i- 
We shall, furthermore, assume that k{t) = k(— r). For data from real detectors the above average 
can be replaced by a time average under the assumption of ergodicity: 



1 /■^/2 

k(t) = - / n{t)n{t - r)dt. (97) 

T J-T/2 



The assumption of stationarity is not strictly valid in the case of real GW detectors; however, 
if their performance doesn't vary greatly over time scales much larger than typical observation 
time scales, stationarity could be used as a working rule. While this may be good enough in the 
case of binary inspiral and coalescence searches, it is a matter of concern for the observation of 
continuous and stochastic GW. In this review, for simplicity, we shall assume that the detector 
noise is stationary. In this case the one-sided noise PSD, defined only at positive frequencies, is 
the Fourier transform of the noise auto-correlation function: 

SnU) ^ ^ / n{T)e^^'f^dT, / > 0, 

J— oo 

= 0, /<0, (98) 

where a factor of 1/2 is included by convention. Since we have assumed that k{t) is an even 
function the above equation immediately implies that Snif) is real. It is quite straightforward 
to show that 

n{f)n*{f') = Sn{mf-f'), (99) 

where n{f) is the Fourier transform of n(t) and n*{f) denotes the complex conjugate of n[f). 
The above identity implies that Sn{f) is positive definite. One derives the above identity by 
expressing the Fourier transforms on the left-hand side by their respective time-domain functions, 
i.e., n{f) = j'^^n{t)e^'^^f^dt, and using Eqs. (|D and 



The autocorrelation function k{t) at r = can be expressed as an integral over Sn{f)- 
Indeed, it is easy to see that 



n^it) = 2 / S„if)df. (100) 
Jo 

The above equation justifies the name power spectral density given to S'„(/). It is obvious that 
Snif) has dimensions of time but it is conventional to use the dimensions of Hz"^ since it is 
a quantity defined in the frequency domain. The square-root of Snif) is the noise amplitude, 
\/Snif), and has dimensions of per root Hz. It is often useful to define the dimensionless quantity 
/i^(/) = fSnif), called the effective noise. In GW interferometer literature one also comes across 
the displacement noise or strain noise defined as heif) = ihnif), and the corresponding noise 
spectrum Siif) = i'^Snif), where i is the arm length of the interferometer. The displacement 
noise gives the smallest strain 6£/i in the arms of an interferometer that can be measured at a 
given frequency. 

^Note that in earlier Sections we have used angular brackets to denote the ensemble average; In this Section, 
however, we shall use an overbar to denote ensemble average while angular brackets will be reserved for denoting 
the scalar product of functions. 
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7.3 Source amplitudes vs sensitivity 

One compares the GW amplitudes of astronomical sources with the instrumental sensitivity and 
assesses what sort of sources will be observable in the following way. Firstly, as comparisons are 
almost always made in the frequency-domain it is important to note that the Fourier component 
h{f) = dt /i(f)e^'^*-^* of a deterministic signal h{t) has dimensions of Hz~^ and the quantity 
f\h{f)\, is dimensionless. It is this last quantity that should be compared with hn{f) to deduce 
the strength of a source relative to detector noise. Secondly, it is quite common also to compare 
the amplitude spectrum per logarithmic bin of a source, \^\h{f)\, with the amplitude spectrum 
of noise, ^/ both of which have dimensions of per root Hz. Justification for these compar- 
isons is given in Sec. |8[ Finally, for monochromatic sources, one compares the effective noise in 
a long integration period with the expected "instantaneous" amplitudes in the following way: 
A monotonic wave of frequency /o observed for a time T is simply a narrow line in a frequency 
bin of width A/ = 1/T around /q. The noise in this bin is £'„(/) A/ = S'„(/)/T. Thus the SNR 
p after a period of observation T is 

^" (101) 



VSn{fo)/T' 

One therefore computes this dimensionless noise spectrum for a given duration of observation, 
Sn{f)/T, to assess the detectability of a continuous GW. 

Surely, if the observation time is T then the total energy (that is, the integrated power 
spectrum) of both the signal and noise must increase in proportion to T^. Then how does the 
SNR for a continuous wave improve with the duration of observation? The point is that while 
the signal energy is all concentrated in one bin, the noise is distributed over the entire frequency 
band. As T increases, the frequency resolution improves as 1/T and the number of frequency 
bins increase in proportion to T. Consequently, the noise intensity per frequency bin increases 
only as T. Now, the signal intensity is concentrated in just one bin since the signal is assumed to 
be monochromatic. Therefore, the power SNR increases as T, or, the amplitude SNR increases 
as VT. 

7.4 Noise power spectral density in first interferometers 

As mentioned in the beginning of this Section the performance of a GW detector is characterised 
by the one-sided noise power spectral density (PSD). The noise PSD plays an important role in 
signal analysis. We shall only discuss the PSDs of interferometric gravitational wave detectors 
since the two prime candidate sources discussed in this Review are both detectable conclusively 
only in a broadband detector. Interferometers have a very broad band sensitivity and are, 
therefore, ideal for the detection of these sources. Since our aim is to foresee the first possible 
detections we shall mainly concentrate on the initial interferometers, GEO, LIGO-I, TAMA and 
VIRGO mentioning future ground- and space-based detectors where appropriate. The sensitivity 
of ground based detectors is limited at frequencies less than a Hertz by the time-varying local 
gravitational field caused by a variety of different noise sources, e.g. low frequency seismic 
vibrations, density variation in the atmosphere due to winds, etc. Thus, for data analysis 
purposes, the noise PSD is assumed to be essentially infinite below a certain lower cutoff fs- 
Above this cutoff, i.e. for / > fg, Table ^ lists the noise PSD 5„(/) for various interferometric 
detectors. The effective noise hn{f) = \/ fSnif) expected in these detectors is summarized in 



Fig. 15. The GEO noise curve is somewhat different from others as it uses a signal enhancement 
technique called 'signal recycling' |165| ] . It is also for this reason that its sensitivity is close to 
LIGO and VIRGO though it is only a sixth/fifth in size as compared to larger interferometers. 

For LISA the Table gives the internal instrumental noise only. However, in the frequency 
range 10~^-3 x 10~^ Hz, LISA will be limited in its sensitivity by the background produced 
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Figure 15: This figure shows the amphtude noise spectral density, \/Sn{f), in initial interferom- 
eters. On the same graph we also plot the signal amplitude, \^\h{f)\, of a binary black hole 
inspiral occurring at a distance of 100 Mpc. Each black hole is taken to be of mass equal to 
10M.Q. (See text in the Section g.2.2 for a discussion.) 
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Figure 16: Same as Fig. but for the LISA detector. Note that this figure, in contrast to 
Fig. 13, uses amplitudes per \/Hz . We also plot signals from supermassive black holes. The 
supermassive BH sources are assumed to lie at a red-shift oi z = 1 but LISA can detect these 
sources with a good SNR practically anywhere in the Universe. 
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Table 5: Noise power spectral densities of initial interferometers, 5'„(/). For each detector the 
noise PSD is given in terms of a dimensionless frequency x = f / fo and raises steeply above a 



Detector 


/s/Hz 


/o/Hz 




104« X 5„(x)/Hz-i 


GEO 


40 


150 


[(3.4x)-3u + 34x-^ + 20(1 -x^ + xV2)/(l + x^/2)] 


LIGO-I 


40 


150 


9.00 [(4.49x)-56 + 0.162;-''-52 + 0.52 + 0.32x2] 


TAMA 


75 


400 


7.50 [x-5 + 13x-i + 9(1 + x^)] 


VIRGO 


20 


500 


3.24 


;(6.23x)-5 + 2x-^ + l + x2] 


LISA 


10-5 


10-3 


420 


(x/5.62x)-i^/3 _^ io3 + 



by several populations of Galactic binary systems, such as closed white-dwarf binaries, binaries 
consisting of cataclysmic variables, etc. This binary confusion noise has been well-modeled 
(see |]l6| and Section ^) and gives the following 'external' LISA 

lO-^^Se(x) = (x/5.56 X 10^)-^ ^ X < 10-°-i5 

= (x/12.3)-^-5, 10-°-15 < ^ < 

= (x/471)-2-6, 10°-25<x, (102) 

where x = f / fo and fo = 10~^ Hz. 

7.5 False alarms and detection threshold 

Gravitational wave event rates in initial interferometers are expected to be rather low: about a 
few per year. Therefore, one has to set a high threshold so that the noise generated false alarms 
mimicking an event is negligible. For a detector output sampled at 1 kHz and processed through 
a large number of filters, say 10^, one has ~ 3 x 10^^ instants of noise in a year. If the noise 
is Gaussian then demanding that no more than one false alarm occurs in a year's observation 
sets the threshold to be about 7.5 times the standard deviation of noise. Therefore, a source 
is detectable only if its amplitude is significantly larger than the effective noise amplitude: i.e. 
fh{f) ^ hn{f). The reason for accepting only such high-sigma events is because the event rate 
of a transient source, i.e. a source lasting for a few seconds to mins, such as a binary inspiral, 
could be as low as a few per year and the noise generated false alarms, at low SNRs ~ 3-4, 
over a period of an year, tend to be quite large. Setting higher thresholds for detection helps in 
removing spurious, noise generated events. However, signal enhancement techniques (cf. Sec. ^) 
make it possible to detect a signal of relatively lower amplitude provided there are a large number 
of wave cycles and the shape of the wave is known accurately. 

7.6 Beam pattern functions 

Gravitational wave detectors are sensitive to waves coming from almost any direction in the sky 
although the degree of their sensitivity depends on the actual direction. The sensitivity of a 
detector to the direction of the wave is described by what is called the beam pattern function, 
also referred to as the antenna pattern j^. A GW antenna responds best if the waves are incident 
at right angles to the principal direction, as in the case of a cylindrical bar antenna, or principal 
plane, as in the case of an interferometric detector. For waves incident in any other direction the 
response will alter by trigonometric factors that will be different for the two polarisations, /i+ 
and /ix. Denoting the two beam pattern functions by -F+ and -Fx, the response C of an antenna 
to a GW of polarisation amplitudes /i+ and /ix is given by 



C 



Fl{e,^,iP)hl + F^{e,ip,i:)hl ^'^ . (103) 
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where {0, if) denotes the dh'ection to the source and if) is the polarisation angle. For a resonant 
bar detector, with its longitudinal axis aligned along the z-axis, the response is, 



F+ = sin^ e cos 2^, Fx = sin^ 6 sin 2^1, (104) 

and for an interferometer with its arms in the (x, y) plane and at right angles to each other and 
the X-axis bisecting the two arms 166 1 



F^{6,ip,il}) = — (1 + cos^ ^) cos cos 2?/^ — cos ^ sin sin 2'i/', 

FxiOjifjip) = -(1 + cos^ 6*) cos 2(^ sin 2-0 + cos sin 2(^ cos 2V'. (105) 

If a gravitational wave lasts long enough, then the detector's motion relative to the source 
will cause the detector to see different polarisations at different times (i.e. C is a function of 
time). Thus, long duration observation will help in resolving the two polarisations. However, 
this will not be possible for inspirals observed by ground-based antennas as the waves will only 
last for a few mins, during which C essentially remains a constant. In this case an interferometric 
antenna will observe only a certain combination of the two polarisations and information of the 
source direction and wave polarisation can only be extracted if several widely separated antennas 
observe the same signal. 



8 Data Analysis 

Observing GWs requires a data analysis strategy which is in many ways different from conven- 
tional astronomical data analysis. There are several reasons why this is so: 

• Gravitational wave antennas are essentially omni-directional with their response better 
than 50% of the average over 75% of the sky. Hence our data analysis systems will have 
to carry out all-sky searches for sources. 

• Interferometers are typically broad-band covering 3 to 4 orders of magnitude in frequency. 
While this is obviously to our advantage as it helps to track sources whose frequency may 
change rapidly, it calls for searches to be carried over a wide-band of frequencies. 

• In contrast to EM radiation, most astrophysical GWs are tracked in phase and the SNR is 
built up by coherent superposition of many wave cycles emitted by a source. Consequently, 
the SNR is proportional to the amplitude and only falls off, with the distance to the 
source r, as 1/r. Therefore, the number of sources of a limiting SNR increases as for 
a homogeneous distribution of sources in a flat Universe, as opposed to EM sources that 
increase only as r^/^. 

• Finally, GW antennas acquire data continuously for many years at the rate of several 
mega-bytes per second. It is expected that about a hundredth of this data will have to 
pass through our search analysis systems. Unless on-line processing can be done we cannot 
hope to make our searches. This places huge demands on the speed of our data analysis 
hardware and a careful study of our search algorithms with a view to making them as 
optimal (maximum SNR) and efficient (least search times) as one possibly can. 

Let us first clarify our notation in this Section and recall two important theorems. We shall 
use x{t) to denote the detector output which is assumed to consist of a background noise n{t) 
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and a useful gravitational wave signal h{t). The Fourier transform of a function x{t) will be 
denoted x{f) and is defined as 

/oo 
x{t)e^^'f'dt. (106) 
-oo 

With this definition the inverse Fourier transform is x{t) = J^^x{f)e~^'^^^^^df. The Fourier 
transform of a real function x{t) obeys x{—f) = x*{f). A shift in the time-domain simply 
appears as a constant phase shift in the Fourier domain, i.e. if x{f) is the Fourier transform of 
x{t) then the Fourier transform of x{t — to) is e^'^^^*°-x{f). 

In this Section we use a system of units in which G = c = 1. Thus, for instance, IM© = 
4.925 X 10~^ s and 1 Mpc= 1.08 x 10^^ s. We would like to caution the reader that h{f) is a 
dimensionless gravitational wave amplitude, while in this Section h{f) is the Fourier transform 
of h{t) and has physical dimensions of time. 



8.1 Matched filtering and Optimal Signal-to-Noise Ratio 

Matched filtering is a data analysis technique that efficiently searches for a signal of known 
shape buried in noisy data |167]. The technique consists in correlating the output of a detector 
with a waveform, variously known as a template or a filter. Given a signal h{t) buried in noise 
n(t), the task is to find an 'optimal' template q{t) that would produce, on the average, the best 
possible SNR. In this review we shall treat the problem of matched filtering as an operational 
exercise. However, this intuitive picture has a solid basis in the theory of hypothesis testing. The 
interested reader may consult any standard text book on signal analysis, for example Helstrom 
1 167], for details. 



8.1.1 Optimal filter 

Let x{t) denote the detector output. If no signal is present then x{t) is just a realisation of noise 
n(t), i.e. x{t) = n{t), while in the presence of a deterministic signal h(t) it takes the form, 

x{t) = h{t-ta)+n{t), (107) 

where ta is a signal that simply shifts the signal relative to the origin of time. In the case of 
signals h{t) whose frequency changes with time different values of ta correspond to different 
frequencies at t = 0. Since a wave detector can only observe signals in a certain frequency 
band ta could mark the time at which the signal enters the detector's sensitivity band. For 
this reason ta is called the time- of- arrival. Given the time-of-arrival the time at which a signal 
leaves the detector is determined by the parameters characterising the signal, such as the masses 
of the component stars in the case of a binary inspiral signal. Time-of-arrival is an important 
parameter in data analysis as it will usually be unknown and must be determined by observation. 
Measuring arrival times in a network of detectors also helps in determining the direction to a 
source. 

The correlation c of a template q{t) with the detector output is defined as 

/oo 
x{t)q{t)dt. (108) 
-oo 

The purpose of the above correlation integral is to improve the visibility of the signal. The 
following analysis reveals how this is achieved wherein we shall work out the optimal filter q{t) 
that maximises the statistical average of the correlation c when a signal h[t) is present in the 
detector output. To do this let us first write the correlation integral in the Fourier domain by 
substituting for x{t) and q{t), in the above integral, their Fourier transforms x(/) and g(/), i.e.. 
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x(t) = x(/) exp (— 27ri/t) d/ and q{t) = J^^q{f) exp {—27rift) df, respectively. After some 
straightforward algebra one obtains 

where q*{f) denotes the complex conjugate of q{f). In general, c consists of a sum of a two 
terms, a filtered signal S and filtered noise A^: 

c = S + N, (110) 

where S = J^o h{f)q*{f)e^^'^''^df and N ^ J^, n{f)q*{f)df. 

Since n is a real random process, c is also a real random process. If n is specified by a Gaussian 
random process with zero mean then c will also be described by a Gaussian distribution function, 
although its mean and variance will, in general, differ from those of n. The mean value of c is, 
clearly, S - the correlation of the template q with the signal /i, since the mean value of n is zero: 

/oo 
h{f)r{f)e'^'^'^df. (Ill) 
-oo 



The variance of c, that is (c — c)^, turns out to be. 



oo 

2 



(c-c)2 = iV2= Sn{f)\q{f)rdf, (112) 



where Sn{f) is the power noise spectral density defined in Eq. (paj) . Now the power SNR is 
defined as = S'^/N'^ and the amplitude SNR is p. The form of integrals in Eqs. ( |111| ) and 
( |112| ) motivates the definition of the scalar product of functions, which could either be templates 
or waveforms. Given two functions a{t) and b{t) we define their scalar product (a, b) to be 

(«' b)^2 r -f- fa(/)6*(/) + a*(/)6(/)l . (113) 



Snif) 

Recall that Sn{f) is real and positive definite, [cf. Eq. (p9|)], consequently, the above scalar 
product defines a positive definite norm: The norm of a, denoted ||a||, is given by 



\a{f)\' 



Snif) 



1/2 



(114) 



Using the reality of the time-domain function h{t) we can write down the SNR in terms of 
the above scalar product as 

^ (he'^^'f'^Snqf 

P = ^ ^T^- (115) 

{Snq,Snq) 

Now, the scalar product of two functions (a, b) acquires its maximum value when a = b. Applying 
this to the above equation one finds that the template q that maximises p, called the optimal 
template, denoted (joptif), is simply 

gopt(/)=7 ^-y^(^) , (116) 



where 7 is an arbitrary constant. The inverse Fourier transform of Eq.( 116 ) gives the optimal 
template qoptit) in the time-domain. One can see that q'opt(0 is the convolution of the time- 
translated signal h{t) with the inverse Fourier transform of 1/S'„(/). Note that to achieve the 
maximum of the SNR the optimal template has to not only match the shape of the signal but 
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also its time-of- arrival ta (cf: the factor e^'^*-^*" in the expression for the optimal template). Since 
one would not know the time-of-arrival of the signal before hand one will have to construct the 
correlation of the detector output for several different relative lags of the template with respect 
to the detector output. In other words, one constructs the correlation function c(t'), 

x{t)q{t-t')dt= / x(/)r(/)e"'"*^* 4f (117) 

-oo J —oo 

where t' is called the lag parameter. 



8.1.2 Optimal signal-to-noise ratio 



We can now work out the maximum, or optimal, SNR by substituting Eq. ( 116 ) for the optimal 
template in Eq. (|115|), 



/Copt = {h, h) 



1/2 



df 



Hf) 



Snif) 



1/2 



(118) 



We note that the optimal SNR is not proportional to the total energy of the signal, which is 
4 Jq°° ci/|/i(/)p, but rather the integrated signal power weighted down by the noise PSD. This 
is in accordance with what we would guess intuitively: the contribution to the SNR from a 
frequency bin where the noise PSD is high, and hence less reliable, should be smaller than from 
a bin where the noise PSD is low. Thus, an optimal filter automatically takes into account the 
nature of the noise PSD. In this final expression for the optimal SNR the parameter ta does not 
appear because the time-of-arrival optimal template matches with that of the signal and hence 
cancels out in the scalar product. 



The expression for the optimal SNR Eq. (118) suggests how one may compare signal strengths 
with the noise performance of a detector. Note that one cannot directly compare |/i(/)p with 
Snif) ^ they have different physical dimensions. In GW literature one writes the optimal SNR 
in one of the following equivalent ways 



Popt 



df 



VfHf) 



f Snif) 



1/2 



df 



fHf) 



f fSnif) 



1/2 



(119) 



which facilitates the comparison of signal strengths with noise performance. One can com- 
pare the dimensionless quantities /|/i(/)| and \/ fSnif) or dimensionful quantities y/f\hif)\ and 

VSM)- 

Signals of interest to us are characterised by several (a priori unknown) parameters, such as 
the masses of the component stars in a binary, their intrinsic spins, etc., and an optimal filter 
must agree with both the signal shape and its parameters. A filter whose parameters are slightly 
mis-matched with that of a signal can greatly degrade the SNR. For example, even a mis-match 
of one cycle in 10^ cycles can degrade the SNR by a factor 2. 

When the parameters of a filter and its shape are precisely matched with that of a signal 
what is the improvement brought about as opposed to the case when no knowledge of the signal 
is available? Matched filtering helps in enhancing the SNR in proportion to the square-root of 
the number of signal cycles covered in the detector band, as opposed to the case when the signal 
shape is not known and all that can be done is to Fourier transform the detector output and 



compare the signal energy in a frequency bin to noise energy in that bin (see, e.g., |168| for a 
proof). 
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8.1.3 Matched filtering of continuous GW 



We will now apply the matched filtering theorem to observations of sources emitting continu- 
ous gravitational waves (CW) at a single frequency. If we observe a monochromatic source of 
frequency /o, that is h{t) = /iq cos(27r/ot), for a duration T then its Fourier transform is a sine 
function (sincx = s'mx/x): 



= (W2)sinc[27r(/-/o)r]. 



(120) 



In the above expression we have ignored the sine function that occurs at / = — /o, as it would 
not contribute to the SNR since the integral in Eq.( |ll8| ) runs only over positive frequencies. 
Since the sine funciton is strongly peaked around f = fo and also because the noise PSD is 
slowly varying in the frequency range 100-1000 Hz, one can treat the power spectrum appearing 



in the SNR integrand Eq.(118) as a constant and write the optimal SNR as: 



CW 

Popt 



V Snifo 



\Hf)\'df 



1/2 



hn 



VSnifo)/T 



(121) 



where we have used the identity /q°° sine {2t: fT)df = 1/AT. Equation (|121|) justifies Eq.(lOl) 
which was derived heuristically. 



8.2 Matched Filtering Inspiral Waves from Compact Binaries 

As pointed in earlier Sees., the last few mins in the evolution of a compact binary (NS-I-NS, 
NS-I-BH, BH-I-BH) is a promising source for interferometers that are currently being built. The 
waveform from these sources is known very accurately and therefore matched filtering is the 
best choice for detecting these sources. Matched filtering is very sensitive to the phasing of the 
waves. It is important, therefore, to keep accurate phase information in our search templates; 
their amplitudes can be taken to be that given by the lowest order post-Newtonian theory. Such 
an approximation which works only with phase corrections and neglects amplitude corrections 



is called the restricted post- Newtonian approximation [16£] 



8.2.1 Accurate templates for inspiral search 

To compute the waveform one must know the evolution of its phase and amplitude, which 
involves the computation of the relative velocity va [t) of the two stars and the phase evolution 
^A{t) of the binary. The subscript A (which will be dropped, for convenience, from va) denotes 
the fact that we know these quantities only in a certain approximation. In the restricted post- 
Newtonian approximation the binary inspiral waveform is given by (see, e.g., Ref. |170f| ) 

h^{t) = hQv\{t) cos ^A[vA{t)\. (122) 

The amplitude Hq depends on the masses Mi and M2 of the component stars, the distance to 
the binary r and the orientation of the source relative to the antenna. More precisely, 

ho = ^C{i,e,^,ij). (123) 
r 

Here, M = Mi -|- M2 is the total mass of the binary, rj = MiM2/M^ is a dimensionless 
(symmetric) mass ratio, which takes a maximum value of 1/4 when Mi = M2. can also be 
thought of as a measure of how strongly the geometry is different from Schwarzschild geometry 
of a single body and is sometimes refereed to as the deformation parameter. C is a function 
of the various angles describing the polarisation of the wave and the orientation of the source 
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relative to the antenna. For a binary whose orbit is inchned at an angle i to the line of sight, C 
is given by 



C{i,e,ip,i;) = VA^ + B^, A= -{l + cos'^i)F+{e,ip,iP), B = cosi F^{e,ip,xP). (124) 

In the above equation the angles 9, ip and ip parameterize both the propagation direction and 
the polarisation of the GW with respect to the detector and and Fx are detector beam 
pattern functions given in Eq.(105). 

For the purpose of computing the SNR from candidate binaries we can either assume the 
source to be ideally oriented, that is 



= V = 0, F+ = l, Fx = 0, C = l, 



(125) 



giving us the best possible SNR /Oideah or we can aim at computing the rms (root-mean-square) 
SNR prms by averaging over all angles. In the latter case we have following rms values 



Fi 



1/2 



F^ 



1/2 / 2^/2 ^2^ 



(126) 



which will be used in the next Section to compute the various SNR. 

The relative velocity VA{t) and the GW phase 4)A{t) = 2ipA{t) in Eq.( |122| ) are given in terms 
of the following ordinary differential equations yi71| : 



:FAiv) 

dt ' ME'^{v) 



dv 

-- + 



0, 



d(t)A 
dt 



2v^ 



0, 



(127) 



V = 



where Ea{v) is the relativistic total energy per unit mass, i.e., E = (Ftotai — M)/M, E'^ 
dEA/dv being its ?;-derivative, and J-'a = —dEA/dt is the GW flux escaping the system, at the 
given approximation]^. The above differential equations can be formally solved to obtain 



tA = to-M 



dv ^ 



v) 



auv'^ 



(128) 



In general, there is no closed-form solution to the above integrals. However, while working 
with post-Newtonian theory one essentially has Taylor series expansions of GW energy and 
flux functions and the above integrals, therefore, can be solved to obtain a solution of the 
form tA = tA{v) and (pA = (pAi^), which can in turn be inverted to obtain va = VA{t) and 
4'A = (pAit)- However, there is much debate in the literature on what is the best representation 
of an exact inspiral wave emitted by a binary system: Some authors work with post-Newtonian 



expansions of flux and energy and the closed form solution mentioned above |172], others with 
the ODEs Eq. (127), or, equivalently, with the integrals Eq. (128) [173, 171], and yet others 
with a more accurate representation of the energy and flux functions, called P-approximants 
| 17C , |174| . These latter representations have been shown to be extremely accurate in the test 
mass limit (i.e. r] ^ 0) and are also expected to be well-behaved when rj ^ |170, 174, 171|. 

At the lowest post-Newtonian order, that is at the quadrupole approximation, and for circular 
orbits we have 

En{v) = -'nv^/2, TNiv) = -32r]'^v^^/5, (129) 



where a subscript is used to denote that the quantity is given at the lowest, i.e. Newtonian, 
order in post-Newtonian theory. Substituting these expressions in Eq.(|128|) and inverting the 



See Al and A_3 for the post-Newtonian expansions of the energy and flux functions. 
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resulting equations gives the following waveform: 

tiv) = to + ^{v^'-v~^) (130) 

m = ^0 + {vo' - v~') . (131) 
h{t) = /iow^(t)cos(/)(t). (132) 

where we have dropped subscript for brevity, vq and (po are the relative velocity of the two 
stars and the phase of the wave at time to and v in the last two equations is that which is 
obtained by inverting the first equation. 

8.2.2 Signal-to-noise ratios for binary inspiral signals 



The optimal SNR Eq. ( |118| ) is given in terms of the Fourier transform of the signal. For compact 
binary inspiral events an estimate of the optimal SNR can be made by using the stationary phase 
approximation to the Fourier transform of the signal [|175|]. In this approximation, the modulus 



of the Fourier transform of the signal in Eq.( 122 ) is found to be [ I74 | 



5 xl/2 f-7/6 



Substituting for /iq in terms of the quantity C and using the resulting expression in Eq. ([TT 
we have the following optimal SNR 



/lso /-7/3]1/2 



df 



is SnU) 



(134) 



where we have now introduced a specific lower and upper limit in the integral. The lower cutoff 
is dictated by the noise PSD that rises steeply below a certain frequency fs- There is hardly 
any SNR to be gained by extending our templates to frequencies below this lower cutoff. For 
ground-based interferometers fs is in the range 1-50 Hz while for the space-based LISA fg is 
slightly smaller than 0.1 mHz. The upper cutoff /lso is determined by the location of the last 
stable orbit of the binary. In the test mass limit (i.e., r/ — > 0) the inspiral signal would terminate 
when the test mass is at a distance -Rlso = 6M from the central body. This corresponds to a 
GW frequency of /lso given by 

^'so = ^^2^ = ^3^' (135) 

where the first equality is just a statement of Kepler's law and in the second equality we have 
used i?LSO = 6M. This upper limit[| is important since it is this which limits the kind of 
inspiral signals that an interferometer can observe. For instance, a ground-based interferometer 
will not be able to observe the inspiral of massive binary black holes of mass ~ 10^-10^ Mq, 
since this signal would terminate at frequencies in the milli-Hertz region, where ground-based 
interferometers are too noisy. However, a space-based detector, such as LISA, is sensitive to the 
inspiral of massive black holes occurring anywhere in the Universe. 

We are finally in a position to write down the rms and ideal SNRs. For a binary at a distance 
r from Earth consisting of stars of individual masses mi and m2, total mass m and symmetric 

^The justification for using an upper limit in the frequency domain based on an upper limit that occurs in 



the time-domain is somewhat technical; the interested reader is referred to [171] and references therein for a 
discussion. 
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Table 6: Signal-to-noise ratios for some archetypal binaries consisting of neutron stars (NS) of 
1.4 Mq and/or black holes (BH) of 10 and 15 Mq at 100 Mpc and for the hypothetical MACHO 
BH binaries of 0.5Mq at 10 Mpc. 



Detector\ Binary 


(0.5,0.5) 


(1.4,1.4) 


(1.4,10) 


(10,10) 


(15,15) 


GEO 


2.6 


0.58 


1.2 


2.8 


3.6 


LIGO 


3.8 


0.95 


2.0 


4.6 


5.5 


VIRGO 


3.5 


0.78 


1.5 


3.4 


4.3 



mass ratio r/ the rms and ideal SNRs are obtained by using the rms and ideal values of C from 



Eq. (|12|) and Eq. (|25D, respectively, in Eq. (|134| ): 



- ^ [ye 



/lso f-7/3nV2 
df- 



Adeal — -Prms- (136) 



Note that the SNR depends only the combination A4 = mrf'^^ of the two masses called the chirp 



mass [176 1 and on the integral whose value crucially depends on the upper limit. In ground- 
based detectors for /lso larger than about 300 Hz the upper limit is not too important while 
for lower values of the LSO frequency the integrals begins to degrade. 

Though the SNR depends only on the chirp mass it does not mean that we cannot measure 
the two masses by GW observations. Recall that the optimal SNRs are obtained only when the 
phases of the signal and the template are finely matched. In matching the phase one resolves 
the degeneracy in the two masses. Indeed, post-Newtonian theory offers the opportunity to 
measure more than two parameters although, in general relativity, each of those parameters, 
in the point particle limit (i.e., neglecting intrinsic spins) depend only on the two component 
masses. This can be potentially used to test general relativity in the strongly non-linear regime. 



(See Ref. [|177| , |178| for a discussion.) 



A detector's sensitivity is sometimes measured by the distance r^'^'^''^^ to an ideally oriented 
NS-NS source that would produces an SNR of 5. We find that r^^'^'^'^^ = 30 Mpc for GEO, 
45 Mpc for LIGO, 1.8 Mpc for TAMA and 50 Mpc for VIRGO. 

The integral in the expression for the SNR Eq. ( |136| ) depends on the detector in question via 
the noise PSD Sn{f) and can be computed (numerically) using the expression for noise PSDs 



given in Table ^. The SNRs so obtained are plotted in Fig, 17 as a function of the total mass 
of a binary consisting of equal masses (i.e. rj = 1/4) and listed in Table ^ for three archetypal 
binaries at 100 Mpc and hypothetical MACHO binaries at 1 Mpc. We note that an inspiral 
event from a BH-BH binary at 100 Mpc will have a reasonably good SNR in all the three larger 
interferometers, namely GEO, LIGO and VIRGO. A coincidence search for these sources can, 
therefore, unambiguously pick out these signals and hence we believe that the first signals to be 
detected are binary black hole inspirals. 

For intermediate mass black hole binaries of mass in the range of 100 solar masses, VIRGO 
and GEO can achieve excellent SNRs. There is no convincing evidence for the existence of 
BHs of such masses, and their binaries, but future GW observations should shed light on these 
systems. 

As mentioned earlier, enhanced interferometers will have ten times better sensitivity and 
bandwidth and will therefore be able to achieve similar SNRs at 1 Gpc. As we have seen in earlier 
Sections at such distances the event rate builds up to several per week. Construction of these 
advanced interferometers will, therefore, make routine detection of compact binary inspirals and 
coalescences possible and thereby open up the new gravitational window for observation. 

The SNRs achievable for the LISA detector is phenomenal and is listed in Table |^ LISA 
should be able to observe ~ million solar mass BH inspiral almost anywhere in the Universe. 
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Figure 17: Signal-to-noise ratio, in initial interferometers, as a function of total mass, for inspiral 
signals from binaries of equal masses at 100 Mpc and averaged over source inclination and 
location. (TAMA has been left out as the SNRs in that case are too low for these sources.) 

The time spent by a binary in the LISA band width can be larger than the life time of the 
detector. For this reason we have expressed the SNR in the above equation in terms of the 
frequency range (/i,/2) in which the signal is extracted during a given observational period. 
For a given duration of observation the binary signal will be the strongest if the observation 
starts closest to the instant of coalescence and therefore we shall take /2 = /lso • 

Table 7: Signal-to-noise ratios in LISA for supermassive binary black holes and stellar mass 
compact objects falling into them, at a cosmological distance of 3 Gpc. 



mi (Mo) 


msCM©) 


/i (Hz) 




/lso (Hz 


) 


P 


10^ 


10^ 


1.08 X 10" 


-5 


2.2 X 10" 


-4 


1900 


10^ 


10^ 


4.54 X 10" 


-5 


2.2 X 10" 


-3 


4900 


10^ 


10^ 


1.92 X 10" 


-4 


2.2 X 10" 


-2 


1700 


10^ 


10 


2.75 X 10" 


-3 


4.4 X 10" 


^3 


21 


10^ 


10 


5.55 X 10- 


-3 


4.4 X 10" 


-2 


11 


10^ 


1.4 


3.93 X 10" 


-3 


4.4 X 10" 


^3 


3.8 


10^ 


1.4 


1.15 X 10" 


-2 


4.4 X 10" 


^2 


2.2 



8.3 Sensitivity to stochastic gravitational waves 

In Sees. |5| and |6| we have seen how astrophysical processes can generate a stochastic back- 
ground of GW and what physical interactions can produce a background in the early Universe, 
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respectively. A stochastic background is easily detected if the instrumental and environmental 
background noise levels in a detector are much smaller than the GW background and especially 
if the statistical properties of the GW background or its spectral characteristics, are different 
from what is expected for the background noise. 

Even when the GW background is far too small compared to the instrumental noise back- 
grounds, the theory of matched filtering developed in the previous Sections, for the detection of 
deterministic signals, can be generalised to the case of a stochastic background of gravitational 
radiation. As opposed to deterministic signals we will not have the advantage of using a template 
that we can compute before hand. Rather, the idea is to use the (noisy) data in one detector, in 
which a stochastic background may be present, as a template to detect the background in (noisy) 
data from another detector. In other words, we can use the 'noisy template' from one detector 
to match filter the background in another. Because our template is noisy, the enhancement in 
SNR will not be as good as in the case of matched filtering a deterministic signal. However, as 
the background will always be present, if we can correlate the outputs of two detectors for a 
long enough duration, in principle, the background will show up above other noise sources. 

The idea that one detector output can be used as a template to enhance the visibility 
of the background in another essentially assumes that (1) the two detectors record identical 
GW signals and (2) the instrumental noise in the two detectors are uncorrelated. In reality 
neither of these assumptions will hold good perfectly. If a pair of detectors happen to be 
close together on Earth, so that they record identical stochastic signals, then they will also be 
affected by similar environmental disturbances such as seismic activity, wind, storm, etc., so 
that the background noise might have a large non-GW correlation. If the two detectors are 
far apart then their environmental disturbances are unlikely to be correlated but they may be 
registering different polarisations and phases of the background GW so that the cross-correlation 
might be insignificant. Clearly, a compromise is in order: Two widely separated detectors are 
good for cancellation of the noise background, while two nearby detectors are favoured for 
enhancing the signal background. In what follows we will derive the SNR enhancement achieved 
on cross-correlating data from two identical detectors placed at the same location. The following 



derivation is a heuristic one to indicate how the cross-correlation works; see Ref. [15S, 16C] for 
a more rigorous derivation involving two detectors of different orientations. 

Let xi{t) and X2{t) denote the outputs of two interferometric antennas located nearby to 
each other and having the same orientation. They will both have the same stochastic signal h{t) 
and we will assume that the internal noises ni{t) and 71-2 (t) are uncorrelated. We will also assume 
that the rms value of the internal noise is much larger than the rms value of the stochastic signal 
we wish to detect. Thus, we can write 

xi{t) =ni{t) + h{t), X2{t) =n2{t) + h{t), ^ = ^ > 7?, (137) 

where an overline indicates the ensemble average and we have assumed that the ensemble av- 
erages of the background noise and stochastic signal are both zero. In order to extract the 
stochastic background we construct the cross-correlation integral C. The cross-correlation of the 
two detector outputs each lasting a time T is given by 

pT/2 

C^= xi{t)x2{t)dt = {xi,X2), (138) 

J-T/2 

where we have introduced the bracket notation to denote the cross-correlation of two data sets 
xi and X2- The cross-correlation is denoted as since it is quadratic in the useful signal h. 
Using Eq. ( |137D can be written as 

= (ni, na) + (ni, /i) + (h, na) + (h, h) . (139) 
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is a random process whose average, under the assumption that the internal detector noises 
are uncorrelated, gives the useful signal S'^ 

= C^ = Jhji). (140) 

From the above equation we conclude that the signal component of grows in proportion to 
the integration time T, i.e., oc T. This, of course, is only the statistical average of and it 
must be compared with the fluctuation in to deduce its significance. Let us define the noise 
component of by 

N"^ = C"^ -C^ = {ni,n2) . (141) 

A^^ is also a random variable and its rms value will involve the noise spectral densities [cf. 
Eqs. (^) and ( |100| )] Smif) and S'n2(/) of the two detectors. Assuming that the two detectors 
have identical noise spectra, say Sn{f), the rms value of N"^ turns out to be 

roo 

m^T dfSlU). (142) 



We see from Eqs. (140) and (142) that the SNR for the stochastic background grows as S/N oc 
T"^/^. Therefore, in principle, a stochastic background that is below the internal noise can be 
extracted by integrating for a sufficiently long duration. Conversely, given that we have data sets 
of a certain duration T with internal noise spectrum Snif), we can estimate the minimum level of 
the stochastic background we can extract. These estimates show that initial interferometers will 
be able to detect a stochastic background at the level of h{f ) ^ 2 x 10~^^, or Qqw ~ 5 x 10~^, 
(cf. Sec.^, with a 90 % confidence level, after three months of integration |160|. Advanced 
ground-based detectors should be able to reach a level that is five orders of magnitude lower (in 
terms of r^cw)- The space-based LISA will have sensitivity to primordial background similar 
to the advanced ground-based detectors but in a frequency range ~ 10^^-lO^^Hz where the 
primordial background has larger spectral amplitudes h{f). Since LISA is likely to be in place 
around the same time as the third generation advanced ground-based detectors, there is an 
exciting possibility to detect primordial background, of the kind discussed in Sec. within the 
next decade. They will also help us to understand populations of astrophysical sources that 
generate a stochastic GW background of the sort discussed in Sec. ^ 



8.4 Computational Costs 

Matched filtering places stringent demands on the knowledge of the signal's phase evolution 
which depends on two things: (1) our modeling of the signal and (2) the parameters characteris- 
ing the signal. If our signal model is inaccurate or if the signal's parameters are unknown, there 
could be a loss in the SNR extracted. For instance, in the case of inspiral signals, a mismatch 
of one cycle in 10^ cycles leads to a drop in the SNR by more than a factor of two, losing a 
factor of eight in the number of potentially detectable events. (Recall that the SNR is inversely 
proportional to the distance to a source; thus an SNR loss by a factor a will reduce the span of a 
detector by the same factor, resulting in a decrease in the volume of observation, and hence the 
number of events, by a factor a^.) Moreover, since the parameters of a signal will not be known 
in advance, it is necessary to filter the data with a family of templates located at various points 
in the parameter space — e.g., placed on a lattice — such that any signal will lie close enough to 
at least one of the templates to have a good cross-correlation with that template. The number 
of such templates is typically very large. This places a great demand on the computational 
resources needed to make an on-line search, that is to search for signals in the detector output 
at the same rate at which data is acquired. We shall discuss below the method of finding the 
number of templates to filter any known signal and the computational resources required to 
analyse the data on-line. In the next Section we will discuss the tools needed in parameter 
estimation. 
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8.4.1 Ambiguity function 



Ambiguity function, well known in statistical theory of signal detection |167], is a very powerful 
tool in signal analysis: It helps in making estimates of variances and covariances involved in 
the measurement of various parameters, in computing biases introduced in using a family of 
templates whose shape is not the same as that of a family of signals intended to be detected, in 
assessing the number of templates required to span the parameter space of the signal, etc. We 
discuss parameter estimation in this Section and computational costs in the next. 

We begin by defining a normalised waveform. A waveform a is said to be normalised if 

1 /2 

it has unit norm: ||a|| = (a, a) ' = 1. The norm of a signal h is also its optimal SNR: 
Popt = {h, h)^^'^ = \ \h\ \ = /iQ. For this reason the norm of a signal is also referred to as the signal 
strength. 

The ambiguity function is defined as the scalar product of two normalised waveforms max- 
imised over the initial phase of one of the waveforms. More precisely, it is the absolute value of 
the scalar product of two normalised waveform^. 

Let a{t; a) be a normalised waveform. Note that we use the symbol a to denote a family of 
waveforms all having the same functional form but differing from one another in the parameter 
values. Indeed, a = {a^\A = 0, . . . ,p} denotes the parameter vector comprising of p + 1 param- 
eters. It is conventional to choose the parameter to be the lag t', which simply corresponds 
to a coordinate time when an event occurs and is, therefore, called an extrinsic parameter, while 
the rest of the p parameters, a^, k = 1, . . . ,p, are called the intrinsic parameters and charac- 
terise the GW source. Given two normalised waveforms, a{t; ai) and a{t;a2), whose parameter 
vectors are not necessarily the same, the ambiguity A is defined as 

A{ai, 02) = niax | (a(ai), 0(02)) | • (143) 

Since the waveforms are normalised it follows that A{ai, ai) = 1, and, if ai is not equal to a2, 
A{ai,a2) < 1. 

It is important to note that in the definition of the ambiguity function there is no need for 
the functional forms of the template and signal to be the same; the definition holds good for 
any signal-template pair of waveforms. Moreover, the number of template parameters need not 
be identical (and usually aren't) to the number of parameters characterising the signal. For 
instance, a binary can be characterised by a large number of parameters, such as the masses, 
spins, eccentricity of the orbit, etc., while we may take as a model waveform the one involving 
only the masses. In the context of inspiral waves a{t; 02) is the exact general relativistic waveform 
emitted by a binary, whose form we do not know, while the template family is a post-Newtonian, 
or some other, approximation to it, that will be used in detecting the true waveform. Another 
example would be signals emitted by spinning neutron stars, isolated or in binaries, whose time 
evolution is unknown, either because we cannot anticipate all the physical effects that affect 
their spin, or because the parameter space is so large that we cannot possibly take into account 
all of them in a realistic search. Of course, in such cases we cannot compute the ambiguity 
function since one of the arguments to the ambiguity function would be unknown. These are 
indeed issues where substantial work is called for. For this Review it suffices to assume that the 
signal and template wave forms are of identical shape and the number of parameters in the two 
cases is the same. 

In the definition of the ambiguity function Eq. ( |143| ) ai can be thought of as the parameters 
of a (normalised) template, while 02 those of a signal. With the template parameters ai fixed, 

^Working with analytic signals h(t) = a(t)e^^'''^*^'^'''°\ where a{t) and </>(f) are the time- varying amplitude and 
phase of the signal, respectively, we observe that the initial phase (j>o of the signal simply factors out as a constant 
phase in the Fourier domain and we can maximise over this initial phase by taking the absolute value of the scalar 
product of a template with a signal. 
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the ambiguity function is a function of p signal parameters a^, k = 1, . . . ,p, giving the SNR 
obtained by the template for different signals. However, one can equally well interpret the 
ambiguity function as the SNR obtained for a given signal by filters of different parameter 
values. Now, the region in the signal parameter space for which a template obtains SNRs larger 
than a chosen value, called the minimal match |17S| ], is the span of that template. Templates 
should be chosen so that together they span the entire signal parameter space of interest with 
the least overlap of one other's spans. 



8.4.2 Metric on the space of waveforms 

The computational costs of our searches and the estimation of parameters of a signal, afford a 
lucid geometrical picture |18[1| , 179|. To develop this picture we begin with the space of signal 



waveforms. A waveform e(t; a), with a given set of values of its parameters, can be thought of 
as a unit vector. (In much of the discussion below we will deal with unit signal and template 
vectors.) The set of sample values h = {hk\k = 0, . . . , — 1} of a waveform, sampled at times 
tk, k = 0, . . . , N — 1, can be thought of as an iV-dimensional signal vector but not necessarily of 
unit norm. Indeed, the output of a detector sampled in the same way can also be regarded as an 
A^-dimensional vector. The set of all detector output vectors forms a vector space. Returning to 
the signal vector, as its parameters are varied the signal vector spans a space in the underlying 
N-dimensional vector space of detector outputs. The dimension of this sub-space is equal to the 
number of parameters of the signal and is called the signal space. 

The signal space, which is a sub-space of the full vector space, has a manifold structure, the 
parameters of the wave constituting a coordinate system and the dimensionality of the space 
being equal to the number of parameters. Having defined a manifold we can ask if it is possible 
to define a meaningful metric on this manifold. Indeed, we already have the necessary tool to 
define the metric, namely the ambiguity function. 

Let us fix the template parameters ai of the template in which case the ambiguity function is 
a function of p+1 signal parameters a^, A; = 0, . . . ,p. Expanding A{ai,a2) about its maximum 
a2 = ai and retaining only quadratic terms we get: 

A{ai,a2) = 1 - 7ABAa^Aa^ + 0[{Aaf], (144) 

where Aq = 02 — ai and 'jab is the metric on the space of waveforms: 

1 d'^A{ai,a2 



7AB(ai) 



(145) 

02=CH 



This is the metric at the point ai on the manifold. From now on we shall drop the suffix 1 on 
the parameter a since it really represents an arbitrary point on the manifold where a template 
resides. Since one can easily maximise over the lag parameter it is desirable to work with 
the metric gi jthat is projected orthogonal to a^, namely 

9ij ^ lij - ^,j = l,---,P- (146) 

Let us suppose that it is required to make a choice of a template bank. By a template bank 
we mean a discrete family of signal waveforms chosen in a given region of the signal parameter 
space. In such a template bank no template waveform will perfectly match an incoming signal 
but if the density of templates is large enough then it may be possible to extract every signal with 
an SNR larger than a certain fraction of the optimal SNR. Most of the early GW data analysis 
literature concerned itself in studying efficient algorithms to set up a template bank so as to 
minimise the computational costs of a search. The density of templates in the bank depends on 
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what is the largest fraction of the optimal SNR one is prepared to lose in a search: Smaller this 
fraction greater is the density of templates. One is normally interested in setting up a bank such 
that each possible signal will have its maximised overlap larger than a certain minimal match 
(MM) with at least some member of the bank. By overlap we mean the scalar product of a 
(normalised) signal waveform with a (normalised) template waveform and maximisation is over 
all the template parameters. The maximised overlap, sometimes referred to as the match, is 
always smaller than or equal to one. Demanding that the proper distance between templates, 
namely gijAa'''Aa^ be as large as possible, for a given minimal match MM, we can obtain the 
following formula for a spacing of templates using Eq. (|144D: 



A,-. = M^^Mm (147) 

V 9kk 

where a factor 2 in the numerator arises because a proper distance of (1 — MM) between the 
least matched signal and a template implies that the proper distance between templates be twice 
that value, assuming that the templates are placed on a square lattice. 

The distance between templates can also be computed more precisely by employing numerical 



methods, as was done in |175| ] (also see [181]). A comparison of the metric-based method 



discussed above and the numerical method shows that the quadratic approximation (147) is 
good typically for A > 0.95. Thus, in the limit of close template spacing, Eq. (|147|) can be used 
to make a choice of templates. 

We now turn our attention to the use of the metric in calculating the number of templates 
needed for a search. If the number M of templates needed to cover a region of interest is 
large, M is well approximated by dividing the proper volume of the region of interest on the 
signal manifold space by the proper volume per template. The proper volume per template, 
AV, depends on the packing algorithm used, which in turn depends on p the dimension of the 
parameter space. For instance, we used a square lattice above. For p = 2, the optimal packing 
is a hexagonal lattice, and thus 

Ay = ^(l-MM). (148) 

There is no packing scheme which is optimal for all p, but it is always possible (though inefficient) 
to use a hypercubic lattice, for which 



AV = {2^{1- MM)/p)P, (149) 

as a starting point. Once the span of a template is known the total number of templates is 
straightforward to compute using 



jd^.Vd.t|i,,ii 

AV ^ ^ 

where p is, as before, the dimension of the parameter space. 

8.4.3 Computational costs for binary inspiral search 

The number of templates for binary inspiral searches has been computed for the post-Newtonian 



signals and ground-based interferometers discussed earlier |181]. The numbers required at a high 
minimal match MM = 0.97 are several hundreds of thousands (cf. Table ^) and imply huge 
computational costs. For instance, to search for inspiral waves on-line (i.e. to search at the 
same rate as the rate at which data is recorded) we would require a dedicated computer that 
can carry out 10^*^ floating point operations per second (i.e., 10 GFLOPS). Building economical 
computers of such speed is the prime concern of various data analysis groups at the moment. 
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Table 8: Number of templates required to search for inspiral waves from a binary consisting of 
stars of mass mi, 7712 > 0.2Mq and minimal match of 0.97. The numbers are given for different 
post-Newtonian families of waveforms. 



Interferometer 


IPN 


1.5PN 


2PN 


LIGO-I 


2.5 X 10^ 


5.3 X 10^ 


4.7 X 10^ 


VIRGO 


1.4 X 10"^ 


1.4 X lO'^ 


1.3 X lO'^ 


GEO600 


4.3 X 10^ 


8.5 X 10^ 


7.5 X 10^ 



8.4.4 Computational costs to search for continuous waves 

The search problem for continuous waves from spinning neutron stars is the most compute- 
intensive job in gravitational wave data analysis. Today, there is little hope that all-sky searches 
lasting for a year or more, can be made. It is easy to see why this is such an intensive job: 
Firstly, the data has to be collected continuously for months together and at a good sensitivity. 
No one has run interferometers for periods as long as that and we do not yet know if this would 
be possible. Secondly, though a neutron star emits a periodic signal in its rest frame, save for 
the neutron star spin-down which indeed induces some modulation in the waveform, because 
of Earth's acceleration relative to the source, the detector does not see a periodic wave. The 
wave is both frequency- and amplitude-modulated. One can, fortunately, de-modulate these 
effects since Earth's motion is known quite accurately, and hence recover the original periodic 
signal. But de-modulation requires a knowledge of the source's direction and its frequency, which 
are unknown in a blind search. The angular resolution one obtains in a year's integration is 
A9 = X/D, where A is the wave length of radiation and D is the baseline of the detector in a 
year's integration, namely 1 A.U. Thus, for / = 100 Hz we have A6 = 10~^ rad or about two 
arcsec. Now, assuming that the source may be in any one of the 4 arcsec^ patches on the sky we 
get the number of patches in the sky for which we will have to try out a de-modulation correction 
to be 47r/(A0)^ = 47rl0^'^. It is quite an impossible task to apply Doppler de-modulation to the 
detector output for each of these ~ 10^^ patches and compute as many Fourier transforms. 

One, therefore, asks the question given a compute power what is the best possible search 
one can do? Is there any advantage in going from a one-step search to a two- or multi-step 
hierarchical search? What about directional searches? These are some of the problems for 
which we have some answer; but a great deal of work is needed and is currently under progress, 
to improve and optimise search algorithms. In the following we will provide a summary of the 
current status. 

The differential geometric formalism discussed above has been used [|182[| to compute the 
number of days of data that a TFLOPS-class computer can analyse on-line (that is, analyse 
T-hours of data in T hours) and carry out a blind (that is, unknown direction, frequency and 
spin-down rate) search. Unfortunately, the longest data we can integrate on-line, for neutron 
stars with spin frequencies / < 100 Hz and spin-down rates less than 1000 years, is about 18 
days. This yields a SNR lower by a factor of 5 as compared to a year's worth of observing. 
On-line searches for neutron stars with / < 500 Hz (largest observed frequencies of millisecond 
pulsars) and spin-down rates of 40 years (shortest observed spin-down rates), can only be made 
for a data set lasting for a duration of 20 hours or less. If the source's position is known in 
advance, but not its frequency, then one can carry out an on-line search, again with a TFLOPS- 
class computer, for the frequency of the source in a data set that is worth 3 months long. This 
is good news since there are many known pulsars and X-ray binary systems that are potential 
sources of radiation. In addition, the obvious targeted search locations are the centre of the 
Galaxy and globular clusters. 
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There have been efforts [183, 184] to study the effectualness of a two-step hierarchical method 
for a bhnd search. Here the basic idea is to construct Fourier transforms of data sets of duration 
smaller than the period in which Doppler modulations will be important and to stack spectral 
densities obtained in this way and to add them all up. This is an incoherent way of building the 
signal since one adds spectral densities that have no phase information. Therefore, one gains 
in SNR less than what an optimal matched filtering method is able to achieve. However, this 
does not matter since (i) the targeted SNR's are quite high ~ 10 and (ii) candidate events can 
always be followed-up using coherent integration methods. These methods afford an on-line all- 
sky blind search (i.e., a search in which no assumptions are made about the parameters of the 
source) for continuous gravitational waves for a period of 4 months or less using a 20 GFLOPS 
computer. Detector groups are planning to build computers of this kind to aid in their search 
for continuous GW. 

8.5 Covariance Matrix and Parameter Estimation 

After a detection has been made, say because a high SNR has been recorded that could not be 
vetoed out, the next step in data analysis is estimation of parameters characterising the event 
and provide error bounds on the measured values. The first thing to note is that one can never 
be absolutely certain that a signal is present in a data train; one can only give confidence levels 
about its presence which could be close to 100% at high values of the SNR. Confidence level is a 
measure of the probability that the observation of an event, such as a large peak in the correlated 
output, is generated by a gravitational wave signal as opposed to a random, non-gravitational 
wave process. The next thing to note is that howsoever high the SNR may be one can't be 
absolutely certain about the true parameters of the signal: At best one can compute a range of 
values in which the true parameters of the signal are most likely to lie. The width of the range 
depends on the confidence level demanded, being larger for higher confidence levels. 

In our search for a signal in the output of a detector we use a discrete, rather than a 
continuous, family of templates. Each template has a particular set of values of the parameters 
and the templates together span an interesting region of the parameter space. The spacing 
between templates in the parameter space will, in general, be quite small. A common estimate of 
the signal parameters is given by the parameters of the template that obtains the maximum SNR. 
Such an estimate is called the maximum likelihood estimate - so-named because the parameters 
of this template maximise what is called the likelihood ratio |167|. Maximum likelihood estimates 



are not always the minimum uncertainty estimates, as has been particularly demonstrated for 
the case of binary inspiral signals |180|. Bayesian estimates, which take into account any prior 



knowledge that may be available about the distribution of the source parameters, as well as the 
knowledge from the output of a whole bank of templates rather than a single template, often 



give a much better estimate [185|. The reason for the better performance of Bayesian estimates 
is that they make a quantitative use of the information at hand. 

In a measurement process any estimation of parameters, howsoever efficient, robust and 
accurate, is unlikely to give the actual parameters of the signal since, at any finite SNR, the 
presence of noise alters the input signal. In geometric terms, the signal vector is being altered 
by the noise vector resulting in a vector that lies outside the signal manifold. Techniques such 
as matched filtering aim at computing the best projection of this altered vector onto the signal 
space. The true parameters of the signal are expected to lie within an ellipsoid of p dimensions at 
a certain level of confidence - the volume of the ellipsoid increasing with the level of confidence. 
The axes of the ellipsoid are the 1-a uncertainties in the estimation of parameters and the 
confidence level corresponding to a 1-a uncertainty is 0.67^, confidence level corresponding to 
a 2-0" uncertainty is 0.95^, and so on. 

The topic of parameter estimation deserves a much wider discussion than given here. How- 



79 



Table 9: Errors in the estimation of instant of coalescence Atp, phase at coalescence A(/>c, 
chirp mass AM/M, and symmetric mass ratio Arj/r], and percentage bias in the estimation of 
the total mass Bm = 100(1 — rnF jm?^^ and percentage bias in the estimation of the mass ratio 
Brj = WO{l — ri'^/ri^), where T stands for the parameter of the template and X for the parameter 
of exact wave form, for the second-post-Newtonian corrected inspiral waveform neglecting the 
effect of spins. Values are quoted for three archetypal binaries consisting of two 1.4 Mq NSs 
(system NS-NS), a 1.4 Mq NS and a 10 M© black hole (system NS-BH), and two 10 Mq black 
holes (system BH-BH). 



System 


Ate (ms) 


A^c" (rad) 


AM/M (10-^) 


At] /rj 




B^ 


NS-NS 


1.07 


2.94 


0.36 


0.28 


0.214 


0.211 


NS-BH 


1.72 


2.27 


2.20 


0.50 


-6.96 


12.2 


BH-BH 


1.50 


2.19 


5.40 


1.50 


1.40 


0.282 



ever, our goal here is only to assemble the necessary tools from estimation theory for ready use. 



Interested reader can consult the ever-growing literature for further details. (See, Ref. [167| for 



estimation theory and [18C, 186, 187, 188] for applications in GW observations 



8.5.1 Covariance matrix 



The scalar product Eq. ( 113 ) induces the following metric on the signal manifold 180 |: 

' da{a) da{a 



G 



3a* ' da^ 



(151) 



where a{a) is a signal vector of unit norm and q* are the coordinates. The metric Gij is 
defined on the manifold of all the parameters including the lag t' and the initial phase (pQ of 
the waveform. It is easy to show that G projected orthogonal to the initial phase (/>o yields the 
metric 7 in Eq.( |145| ). 

The metric on the space of waveforms introduced above is also called the Fisher information 
matrix |167]. Indeed, it contains the 'information' about how similar or dissimilar are the 
waveforms in a small neighbourhood around the parameter a of the signal. Large values of the 
metric imply that even small changes in signal and template parameters can greatly affect their 
overlap while the opposite is true when the metric coefficients are small. 

The inverse of the information matrix is the covariance matrix Cij, whose diagonal and off- 
diagonal elements are, in the limit of large SNR, the variances in the measured values of the 



parameters and correlation coefficients among different parameters, respectively [167]: 



Gij = (G) 



-1 

ij 



(152) 



Covariance matrix based errors in the estimation of the total mass, reduced mass and the 



instant of coalescence tc have been computed by Poisson & Will ]188] for the second post- 
Newtonian inspiral waveforms. They are listed in Table |^ for three archetypal binaries. The 
relative errors are smaller, for a given SNR, in the case of lighter binaries. This is because lighter 
binaries last for a longer duration and have a larger number of cycles in the detector, making 
it relatively easier, as compared to higher mass binaries, to discriminate waveforms of different 
parameters. 



8.5.2 Biases in estimation 

There are two ways in which an error can occur in the estimation of signal parameters. Firstly, 
an error in the measurement of a parameter occurs because of internal noise which alters the 
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input signal and hence, in the process of maximising the SNR, we err in the estimation of 
parameter by an amount that depends on the SNR. This type of error is a random error and 
normally goes down in inverse proportion to the SNR. In the limit of an infinitely large ensemble 
of measurements the estimated values converge to the true values of the signal parameters. 
Secondly, our estimator may be a biased estimator in the sense that the average over an ensemble 
of measurements might converge to values different from the true values of the parameters. This 
can happen because of one or both of the following two reasons: The search templates we use 
in our detection algorithms might only be an approximation to the true signal, as in the case 
of inspiral wave searches where we use post-Newtonian search templates to look for the fully 
general relativistic signal. As we shall show below, such a search would induce a bias in the 
estimation of parameters, which we shall refer to as the bias of the first kind. Alternatively, 
the estimator might be inherently a biased estimator and may give erroneous values even if the 
search templates are not faulty which we shall refer to as the bias of the second kind. 

An example of the bias of the second kind is the following. In the maximum likelihood 
method we always aim for the largest value of the likelihood ratio. This has the effect of making 
a higher estimation of the amplitude parameter of the signal as shown in [180|. For inspiral 
signals, this would mean that the distance to the source would be under-estimated. 

Biases in parameter estimation for inspiral signals have been discussed in |189|. Biases in the 
estimation of the total mass and the symmetric mass ratio on using standard post-Newtonian 
approximation are quoted in Table |9| and discussed in detail in |170| ] . They find that while using 
standard post-Newtonian approximation to the waveform the bias could be quite large, whereas 
P-approximants to inspiral waves greatly reduce the bias in the estimation of parameters. 



8.6 Conclusions 

Gravitational waves in ground-based interferometers are expected to be below the noise levels of 
first and second generation instruments. We have seen in this Section how a prior knowledge of 
the signal's shape can be used in enhancing the visibility of the signal. With the aid of matched 
filtering a network of initial interferometers, consisting of GEO, LIGO and VIRGO, should be 
able to survey a volume of 3 x 10^ Mpc'^ for binary black holes at a minimum SNR of 3. In 
this volume conservative estimates of coalescence rates predict about a few per year. Thus, the 
binary black holes may be the first events to be registered in our detectors. 

Initial interferometers should also be able detect primordial gravitational wave background 
at the level of Oqw ~ 5 x 10~® by cross-correlating data in nearby interferometer pairs (LIGO- 
LIGO or GEO- VIRGO). However, this is at a level rather too small from a theoretical point 
of view. New data analysis algorithms that exploit the specific signature of the primordial 
background may aid in detecting the background at a much lower level but, at the moment, we 
do not know how to achieve this. 

Initial interferometers will also be able to detect continuous gravitational waves from newly 
born or rapidly spinning non-spherical neutron stars provided that the amplitude is h > 10~^^ 
and the signal lasts for a few months. This, of course, assumes that we know the evolution of 
the phase of the signals emitted in the process which, for old neutron stars is a simple sinusoid. 
However, in the case of waves emitted due to one of the instabilities we do not have the right 
templates. This is also true of transient waves emitted during a supernova explosion or the 
merger of black holes. Thus, much work is needed in understanding sources of gravitational 
waves. Indeed, it may very well be that the new generation of gravitational wave antennas will 
open up a new window for observing the Universe which will aid in our understanding of strong 
and non-linear gravity. 
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1 Appendix 

Here we summarize some basic formulae relevant for the description of massive binary evolution. 



Al Keplerian binary system and radiation back reaction 

Binary stars is one of the main topics of the present paper, so it is necessary to remind the 
reader some basic facts about Keplerian motion in a binary system. The stars are highly 
concentrated objects, so their treatment as point masses is usually adequate for the description 
of their interaction in the binary. Further, the Newtonian gravitation theory is sufficient for this 
purpose as long as the orbital velocities are small in comparison with the speed of light c. The 
systematic change of the orbit caused by the emission of gravitational waves will be considered 
in a separate paragraph below. 



Al.l Keplerian motion 

Let us consider two point masses Mi and M2 orbiting each other under the force of gravity. It 



is well known (see | 190 |) that this problem is equivalent to the problem of a single body with 
mass fi moving in an external gravitational potential. The value of the external potential is 
determined by the total mass M. The total mass M of the system is 

M = M1+M2, (Al) 

and the reduced mass /j. is 

M1M2 , . , 

, = (A2) 

The body moves along an elliptic orbit with eccentricity e and major semiaxis a. The orbital 
period P and orbital frequency O = 2tt/P are related with M and a by the 3d Kepler's law 

This relationship is true for any eccentricity e. 

Individual bodies Mi and M2 move around the barycentre of the system in elliptic orbits 
with the same eccentricity e. Major semiaxes of the two ellipses are inversely proportional to 
the masses: 

and satisfy the relationship a = oi + 02. The position vectors of the bodies from the system's 
barycenter are ri = M2r/ [M1 + M2) and r2 = —Mir/ (Mi +M2), where r = ri — r2 is the relative 
position vector. Therefore, the velocities of the bodies with respect to the system's barycentre 
are related by 

Vi Mo 

and the relative velocity is V = Vi — V2. 

The total conserved energy of the binary system is 

MiVi^ M2V2^ GM1M2 fiV^ GM1M2 GM1M2 , , 

^ = ^— + ^ ^ = -2 ^ = 2^ 
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where r is the distance between the bodies. The orbital angular momentum vector is perpen- 
dicular to the orbital plane and can be written as 



Jorb = MiVi X n + M2V2 X r2 = /iV X r. (A7) 
The absolute value of the orbital angular momentum is 



I Jorbl = /iVGMa(l - e2) . (A8) 

For circular binaries with e = the distance between orbiting bodies does not depend on 
time 

r{t, e = 0) = a 

and is usually referred to as orbital separation. In this case, velocities of the bodies, as well as 
their relative velocity, are also time- independent: 



V = \Y\=na= ^GM/a , (A9) 
and the orbital angular momentum becomes 

|J„.^b| = ^Va = fiQa^. (AlO) 

A1.2 Gravitational radiation from a binary 

The plane of the orbit is determined by the orbital angular momentum vector Jorb- The line 
of sight is defined by a unit vector n. The binary inclination angle i is defined by the relation 
cosi = {n, J orb / Jorb) such that i = 90° corresponds to the system visible edge-on. 

Let us start from two point masses Mi and M2 in a circular orbit. In the quadrupole 
approximation |^], the two polarization amplitudes of GW at a distance r from the source are 
given by 

(-<5/3 1 

/i+ = — J 2(1 + cos i)(7r/M)2/3^cos(27r/t) (All) 

= ±—, 4 cos i(TTfMf/^n sin(27r/t) (A12) 

Here M is the total mass, /x is the reduced mass, and / = rj/vr is frequency of the emitted GW 
(twice the orbital frequency). Note that for a fixed distance r and a given frequency /, the GW 
amplitudes are fully determined by /iM^/^ = Ai^^^, where the combination 

M = ^3/^M2/5 

is called the "chirp mass" of the binary. After averaging over orbital period (so that the squares 
of periodic functions are replaced by 1/2) and orientations of binary orbital plane, one arrives 
at the averaged (characteristic) GW amplitude 

/.(/,Al,r) = ((/^^) + (4))^/^= (f )'^'^^(-/)^/=' (A13) 
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A1.3 Energy and angular momentum loss 



In the approximation and under the choice of coordinates that we are working with, it is sufficient 
to use the Landau-Lifshitz gravitational pseudotensor ||^] when calculating the gravitational 
wave's energy and flux. (This calculation can be justified with the help of a fully satisfactory 
gravitational energy-momentum tensor that can be derived in the field-theoretical formulation 
of general relativity |191| ]). Energy dE carried by a gravitational wave along its direction of 
propagation per area dA per time dt is given by 



dE 



F 



dh+ 
~dt 



2 fdh, 



V dt 



dAdt levrG 

The energy output dE/dt from a localized source in all directions is given by the integral 

dE 



(A14) 



Replacing 



and introducing 



~dt 



dt 



+ 



F{e,(l))r'^dn. 



(A15) 



fdhy 



V dt 
1 



h^{e,(l))dn, 



we write Eq. ( A15| ) in the form 



(A16) 



dE 

m G' 

Specifically for a binary system in a circular orbit, one finds the energy loss from the system 
(sign minus) with the help of Eqs. ( |A16 ) and ( A13| ) : 



dE 
'dt 



^(A^vr/)W3. 



(A17) 



This expression is exactly the same one that can be obtained directly from the quadrupole 
formula [||] 

dE 32 G'^M?M^M , 

-r = -T — r , — (A18) 

rewritten using the definition of the chirp mass and the Kepler's law. Since energy and angular 
momentum are continuously carried away by gravitational radiation, two masses in orbit spiral 
toward each other, thus increasing their orbital frequency 0,. The GW frequency f = 0,/tt and 
the GW amplitude h are also increasing functions of time. The rate of the frequency change is f\ 

Q5/3 



-vr 



8/3_y^5/3y.ll/3^ 



(A19) 



In spectral representation, the flux of energy per unit area per unit frequency interval is 
given by the right-hand-side of the expression 



dE 



dAdf G 2 



i\hifh\' + \Hf). 



G 2 



Slif), 



(A20) 



where we have introduced the spectral density S'^{f) of the gravitational wave field h. In case 
of a binary system, the quantity Sh is calculable from Eqs. (All) and ( [A12| ): 

GV3.^5/3 1 ^^^^^ 



Sh 



c3 12 r2 (7r/)7/3 • 



''a signal with such an increasing frequency is reminiscent of a chirp of a bird. This explains the origin of the 
term "chirp mass" for the parameter A4 which fully determines the GW frequency and amplitude behaviour. 
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A1.4 Binary coalescence time 



A binary system in a circular orbit loses energy according to Eq. (|A17 ). For orbits with non-zero 



eccentricity e, the right-hand-side of this formula should be multiplied by the factor 



(see |192| ]). The initial binary separation oq decreases and, assuming Eq. ( A18| ) is always valid, 
it should vanish in a time 

As we noted above, gravitational radiation from the binary depends on the chirp mass which 
can also be written as = Mrf^^ , where rj is the dimensionless ratio "q = ^x/M. Since r] < 1/4, 
one has M ^ 0.435M. For example, for two NS with equal masses Mi = M2 = 1.4M0 the chirp 
mass is ~ 1.22M0. This explains the choice of normalization in Eq. ( A22| ). 



The coalescence time for an initially eccentric orbit with eo 7^ and separation oq is shorter 
than the coalescence time for a circular orbit with the same initial separation oq [192|: 



tc(eo) = to /(eo) (A23) 



where the correction factor /(eo) is 



eo/i , 121 2 ^ 1181/2299 



19 48/19 ^ 121 2\ 3480/2299 J n_ e2)3/2 

+ 304^0 







To merge in a time interval shorter than the Hubble time {tn ~ 10 Gyr), the binary should have 
a small enough initial orbital period Po ^ -fcr(eO)-A^) and, accordingly, a small enough initial 
semimajor axis < acr(eo,A4). These critical orbital periods and semi-major axes are shown 



as functions of the initial eccentricity eo in Fig. |l^ and 19 , respectively. The lines are plotted 
for three typical sets of masses: two neutron stars with equal masses (1.4M0-I-1.4M0), a black 
hole and a neutron star (IOM0-I-I.4M0), and two black holes with equal masses (IOM0-I-IOM0). 
Note that in order to get a significantly shorter coalescence time, the initial binary eccentricity 
should be eo ^ 0.6. 

A2 Mass transfer modes and mass loss in binary systems 

The gravitational wave emission is the sole factor responsible for the change of orbital parameters 
of a pair of compact (degenerate) stars. However, at the early stages of binary evolution, it is the 
mass transfer between the components and the loss of matter and its orbital momentum that play 
dominant dynamical role. Strictly speaking, these processes should be treated hydrodynamically 
and they require complicated numerical calculations. However, binary evolution can also be 
described semi-qualitatively, using a simplified description in terms of point-like bodies. The 
change of their integrated physical quantities, such as masses, orbital angular momentum, etc. 
governs the evolution of the orbit. This description turns out to be successful in reproducing 



the results of more rigorous numerical calculations (see e.g. [40| for a review). In this approach, 
the key role is allocated to the total orbital angular momentum J^rb of the binary. 

Let star 2 lose matter at a rate M2 < and let /3 (0 < /? < 1) be a fraction of the ejected 
matter which leaves the system (the rest falls on the first star), i.e. Mi = —(1 — (3)M2 > 0. 
Consider circular orbits with orbital angular momentum given by ( |Aiq ). Differentiate both 
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Figure 18: The maximal initial orbital period (in hours) of two point masses which will coalesce 
due to gravitational wave emission in a time interval shorter than lO^'^ years, as a function of 
the initial eccentricity cq- The lines are calculated for IOMq+IOMq (BH+BH) 
(BH+NS), and 1.4Mq+1.4M0 (NS+NS). 
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Figure 19: The maximal initial semimajor axis (in R©) of two point masses which will coalesce 
due to gravitational wave emission in a time interval shorter than 10^*^ years, as a function of 
the initial eccentricity cq. The hues are calculated for IOMq+IOM© (BH+BH), IOM0+I.4M0 
(BH+NS), and I.4M0+I.4M0 (NS+NS). 
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parts of Eq. ( |A10D by time t and exclude dVL/dt with the help of the third Kepler's law (|A3|). 
This gives us the rate of change of the orbital separation: 



a 
a 



l-L. g - I) — ± _ z. — £ 



Mo Jnrh , . 

W + '^T^ (A25) 

M2 Jorb 



One defines the mass transfer as conservative if both /3 = and Jorb = 0. The mass transfer is 
called non-conservative if at least one of these conditions is violated. 

For massive binaries, which we are mostly interested in, it is important to distinguish some 
specific cases (modes) of mass transfer. They amount to: (1) conservative accretion mode, (2) 
non-conservative Jeans's mode (or fast wind mode), (3) non-conservative isotropic re-emission, 
(4) sudden mass loss from one of the components during supernova explosion, and (5) common- 
envelope stage. Separately we consider orbit evolution due to (6) gravitational wave emission, 
which becomes the main factor for short-period compact binaries. For non-conservative modes, 
one can also introduce some subcases, such, for example, as a ring-like mode in which a cir- 
cumbinary ring of expelled matter is being formed (e.g. [193| ). Here, we will not go into the 
details of such subcases. 

A2.1 Conservative accretion. 

In the case of conservative accretion, matter from M2 is fully deposited to Mi. The transfer 
process preserves the total mass M = const (/? = 0) and the orbital angular momentum Jorb = 
const of the system. It follows from Eq. ( [A 2 51 ) that 

MiM2\/a = const , 

so that the initial and final binary separations are related as 

2 

I • (A26) 



af _ MuM2i 



ai \MifM2fJ 

The well-known "rule of thumb" for this case says that the orbit shrinks when the more massive 
component loses matter, and the orbit widens in the opposite situation. During such a mass 
exchange, the orbital separation passes through a minimum, if the masses become equal in course 
of the mass transfer. 

A2.2 The Jeans (fast wind) mode. 

In this mode the ejected matter completely escapes from the system, that is, (5 = 1. The escape 
of matter can take place either in a spherically symmetric way or in the form of bipolar jets 
moving from the system at high velocity. In both cases, matter carries away some amount of 
the total orbital momentum proportional to the orbital angular momentum J2 = {Mi /M) Jorb 
of the mass loosing star (we neglect a possible proper rotation of the star, see jl94| ]). For the 
loss of orbital momentum Jorb it is reasonable to take 

Jorb = ^J2. (A27) 
M2 



In the case /? = 1, Eq. ( A25| ) can be written as 

(Oa^)- Jorb M1M2 



^a? Jorb MM2 ■ ^^^^^ 
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Then Eq. ( A28| ) in conjunction with Eq. ( |A27 ) give J7o^ = const, that is, ^jGaM = const. 

(A29) 



Thus, as a result of such a mass loss, the change in orbital separation is 

aj _ Mi 



m Mf 

Since the total mass decreases, the orbit always widens. 
A2.3 Isotropic re-emission. 

The matter lost by star 2 can first accrete to star f, and then, a fraction j3 of the accreted 
matter, can be expelled from the system. This happens when a massive star transfers matter to 
a compact star on the thermal timescale (< 10^ years). The accretion luminosity may exceed 
the Eddington luminosity limit, and the radiation pressure pushes the infalling matter away 
from the system, in a manner similar to the spectacular example of the SS 433 binary system. 
In this mode of the mass transfer, the binary orbital momentum carried away by the expelled 
matter is determined by the orbital momentum of the accreting star Mi, rather than by the 
orbital momentum of the mass-loosing star M2. The orbital momentum loss can be written as 

Mo 

J orb = Pj^^Ji , (A30) 

where Ji = {M2/M)Jorb is the orbital momentum of the star Mi. In the limiting case when all 
the mass attracted by Mi is fully pushed away, /3 = 1, Eq. ( [A30|) simplifies to 

jpTb ^ M2M2 

Jorb MiM ■ ^ ^ 

After substitution of this formula into Eq. ( A25| ) and integration over time, one arrives at 



The exponential term makes this mode of the mass transfer very sensitive to the components 
mass ratio. If M1/M2 <C 1, the separation a between the stars may decrease so greatly that the 
approximation of point masses becomes invalid. The tidal orbital instability (Darwin instability) 
may set in, and the compact star may start spiraling toward the companion star center (the 
common envelope stage; see section A2.5| below). 



A2.4 Supernova explosion. 

Supernova explosion in a binary system occurs on a timescale much shorter than the orbital 
period, so the loss of mass is practically instantaneous. This case can be treated analytically 



(e.g. 1 195, 196, ^). In general, the loss of matter and radiation is aspherical, so that the 
remnant of the supernova explosion (neutron star or black hole) acquires some recoil velocity 
called kick velocity w. In a binary, kick velocity should be added to the orbital velocity of the 
pre-supernova star. 

The usual treatment proceeds as follows. Let us consider a pre-SN binary with initial masses 
Ml and M2. The stars move in a circular orbit with orbital separation aj and relative velocity 
Vj. The star Mi explodes leaving a compact remnant of mass Mc- The total mass of the binary 
decreases by the amount AM = Mi — Mc- The compact star acquires some kick velocity w. 
Unless the binary is disrupted, it will end up in a new orbit with eccentricity e, major semiaxis 
Of, and the angle 6 between the orbital planes before and after the explosion. In general, the new 
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barycenter will also receive some velocity, but we neglect this motion. The goal is to evaluate 
the parameters aj, e, and 6. 

It is convenient to work in an instantaneous reference frame centered on M2 right at the 
time of explosion. The x-axis is the line from M2 to Mi, the y-axis points in the direction 
of Vj, and the z-axis is perpendicular to the orbital plane. In this frame, the pre-SN relative 
velocity is Vj = (0,1/^,0), where Vi = ^/G{Mi +M2)/ai (Eq. @). The initial total orbital 
momentum is Jj = /ijaj(0, 0, — 1^). The explosion is considered to be instantaneous. Right 
after the explosion, the position vector of the exploded star Mi has not changed: r = (aj,0,0). 
However, other quantities has changed: Vj = {wx-, Vi+Wy, Wz) and J/ = Hfai{0, Wz, —{Vi+Wy)), 



(w 



Xl Wy, 



,Wz) is the kick velocity and jU/ = McM2/{Mc + M2) is the reduced mass 



where w 

of the system after explosion. The parameters o/ and e are being found from equating the 
total energy and the absolute value of orbital momentum at the initial circular orbit to their 
expressions at the resulting elliptical orbit (see Eqs. (A6), ( A10| ), and (|A^)): 



GMcM2 



GMcM2 



2a 



f 



tifai^wl + {Vi + WyY = fif^GiMc + M2)a/(1 - e^) 
For the resulting a/ and e one finds 



(A33) 



(A34) 




Wl+Wl + {Vi + Wy) 



and 



1 2 '^i 
1 - e = X— 



Wl + {V+Wyf 



y.2 



(A35) 



(A36) 



where x = {Mi + M2)/{Mc + M2) > 1. The angle 6 is defined by 



cos I 



which results in 



cos 6 



Vi + Wy 



^Wl + {V^+Wyf 



(A37) 



The condition of disruption of the binary system depends on the absolute value Vf of the 
final velocity, and on the parameter x- The binary disrupts if its total energy defined by the 
left-hand-side of Eq. ([A3^ ) becomes non-negative or, equivalently, if its eccentricity defined 
by Eq. ( A36 ) becomes e > 1. From either of these requirements one derives the condition of 
disruption: 

(A38) 




The system remains bound if the opposite inequality is satisfied. Eq. ( A38| ) can also be written 
in terms of the escape (parabolic) velocity Ve defined by the requirement 



V^ GMcM2 



0. 



Since x = M/{M - AM) and V^ = 2G{M - AM)/ai = 2V^/x, one can write Eq. ( |A38| ) in the 
form 

Vf>Ve. (A39) 



89 



The condition of disruption simplifies in the case of a sphericahy symmetric SN explosion, that 
is, when there is no kick velocity, w = 0, and, therefore, Vj = Vi. In this case, Eq. (^) reads 
X > 2, which is equivalent to AM > M/2. Thus, the system unbinds if more than a half of 
mass of the binary is lost. In other words, the resulting eccentricity 

Ml- Mc 

e = — (A40) 



following from ( A35 ), ( A36| ) , and w = becomes larger than 1, if AM > M/2. So far, we have 



considered an originally circular orbit. If the pre-SN star moves in an originally eccentric orbit, 
the condition of disruption of the system under symmetric explosion reads 

AM = Ml - Mc > -— , 
2 at 

where r is the distance between the components at the moment of explosion. 
A2.5 Common envelope stage. 

This is a very important stage in binary evolution. A possibility of this stage was first suggested 



in |197| . Generally, it occurs in binary systems where the mass transfer from the mass-losing 
star is high, and the companion cannot accrete all the matter. The common envelope stage 
appears unavoidable on the observational grounds. The evidence for a dramatic orbital angular 
momentum decrease at some preceding evolutionary stage follows from observations of certain 
types of close binary stars. They include cataclysmic variables, in which a white dwarf accretes 
matter from a small red dwarf main-sequence companion, planetary nebulae with double cores, 
low-mass X-ray binaries and X-ray transients (neutron stars and black holes accreting matter 
from low-mass main-sequence dwarfs). The radii of progenitors of compact stars in these binaries 
typically should have been 100-1000 solar radii, that is, much larger than the observed binary 
separations. This testifies to some dramatic reduction of the orbital momentum at earlier stages 
of evolution and eventual removal of the common envelope. 

There is no exact criterion for the formation of a common envelope. However, a high mass 
overflow onto a compact star from a normal star is always expected when the normal star leaves 
the main-sequence and develops a convective envelope. The critical mass ratio for the unstable 
Roche lobe overflow depends on specifics of the stars, but is close to 1. Another way for the 
formation of a common envelope is direct penetration of a compact star into the dense outer 
layers of the companion. This can happen as a result of the Darwin tidal orbital instability 
in binaries [|198| , |199(| , or when a compact remnant of supernova explosion with appropriately 
directed kick velocity finds itself in an elliptic orbit whose minimum periastron distance aj(l — e) 
is smaller than the stellar radius of the companion. 



A simplified treatment of the common envelope stage is usually done as follows |200|. The 
orbital evolution of the compact star Mc inside the envelope of the normal star Mi is driven by 
the dynamical friction drag. This leads to a gradual spiral-in process of the compact star. The 
released orbital energy AEorb, or a fraction of it, can become numerically equal to the binding 
energy Ebind of the envelope with the rest of the binary system. It is usually assumed that 
this equality provides a condition for expulsion of the common envelope. What remains of the 
normal star Mi is its stellar core Msc- The final orbital parameters are being derived from the 
condition -Efcmd = acE^Eort, where acE is the efficiency parameters less or equal to one. This 
condition reads 

GMi{Mi-Msc) (GMMsc GMiMA 
XR, = [-2^j 2^ ) ' 
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where aj and a/ are the initial and the final orbital separations, A is a numerical coefficient of 
order 1. Rl is the Roche lobe radius of the normal star approximated by | 201 ] 



Rl 

ai 



0.49 



0.6 + ^2/3 ln(l +g-i/3) 
and q = M1/M2. From Eq. (|A41| ) one derives 



(A42) 



ai Ml 



XoceRl Mc 



Ml AM ' 



(A43) 



where AM = Mi — Mgc- Recent studies (e.g. [p02| , |68| ) show that ac^A falls within the range 
from 0.5 to 2. The mass Msc of a helium core of a massive star is (see j203|] ) 



Mhc ^ 0.073(Mi/M(r 



0J 



a.42 



(A44) 



so the orbital separation during the common envelope stage may decrease as much as by factor 
30-60. 



A3 Post-Newtonian expansions of GW flux and energy 

The gravitational wave flux escaping a system of two compact stars in quasi-circular orbit when 
their orbital frequency is /orb (the dominant GW frequency being / = 2/orb)i is given by 



[ 20^ , |205| , |20q , |207|] 



32rj^v^^ /1247 35r?\ o . 3 

— 1 \ + Airv-^ 

5 |_ V 336 12 J 

/447II 9271r/ 65r/2\ ^ /8191 535?7\ 5' 
I 9072 ^ 504 ^ ) ^ ~ v"672" ^ j 



(A45) 



where v = (vrM/)^/^ is the relative velocity of the two stars and rj = MiM2/M'^ is the symmetric 
mass ratio. The r/-parameter takes a maximum values of 1/4 when the two masses are equal. It 
characterises the extent to which two-body effects in the system are important. The relativistic 
energy E{y ) of the system is given by 



E{v) 



rjv 



9 + ^\.,2 /81-57r/ + 7?2\ ^ 



12 



24 



(A46) 



Solving the differential equations in Eq. ( |127] ) for time and phase in terms of v gives the 
following equations: 



+ 



5M 

/ 3058673 5429r/ 617r/M , /7729 



\ 508032 



+ 



327r 

, H I V" 1 

V252 3 J 5 

,4 
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(A47) 
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/3715 5577 \ 2 3 

V1008 12 y 



( 15293365 ^ 271457/ ^ 3085r/2^ ^4 



\ 1016064 



1008 
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(A48) 
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One can invert the first of the equations above to express v in terms of a post-Newtonian 
expansion in t and then use the resulting expression in the second equation to arrive at an expUcit 
phasing formula. Introducing a new time parameter 6 defined hy 9 = [ri{t]so — t)/{5M)]~^^^ , 
where tiso is a reference time taken to be the time at which the GW frequency is equal to twice 
the orbital frequency at the last stable circular orbit, we find 




(A49) 
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